{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ee8ca02e",
   "metadata": {},
   "source": [
    "# tensorflow\n",
    "\n",
    "这里的代码是用 TensorFlow2.x 完成的，但是应该也可以使用 TensorFlow1.x /PyTorch/MXNet/PaddlePaddle/CNTK 等其他深度学习框架实现。\n",
    "\n",
    "https://www.tensorflow.org/api_docs/python/tf\n",
    "\n",
    "https://www.tensorflow.org/guide/autodiff\n",
    "\n",
    "https://www.tensorflow.org/guide/advanced_autodiff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b4c691bc",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'2.6.0'"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import tensorflow as tf\n",
    "\n",
    "tf.__version__"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "04786b8f",
   "metadata": {},
   "source": [
    "TensorFlow 的 Eager Execution 是一种命令式编程环境，可立即评估运算，无需构建计算图：运算会返回具体的值，而非构建供稍后运行的计算图。这样能使您轻松入门 TensorFlow 并调试模型，同时也减少了样板代码。\n",
    "\n",
    "在 Tensorflow 2.0 中，默认启用 Eager Execution。\n",
    "```python\n",
    "import tensorflow as tf\n",
    "tf.executing_eagerly()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc8ea60f",
   "metadata": {},
   "source": [
    "## 张量操作\n",
    "tf.Tensor\n",
    "\n",
    "https://www.tensorflow.org/api_docs/python/tf/Tensor\n",
    "\n",
    "请百度谷歌：TensorFlow 张量操作"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5099d79e",
   "metadata": {},
   "source": [
    "### 数学操作符 \n",
    "\n",
    "https://www.tensorflow.org/api_docs/python/tf/math\n",
    "\n",
    "不同版本的TensorFlow的算术操作命令并不一样，这很令人恼火。\n",
    "\n",
    "#### TensorFlow的算术操作如下：\n",
    "\n",
    "如果 AttributeError ，尝试  tf.math.xxx，或者 搜索 tensorflow math-ops 或者 tf.math\n",
    "\n",
    "| 操作  [如果 AttributeError ，尝试  tf.math.xxx]                      | 描述                                                         |\n",
    "| :-------------------------- | :----------------------------------------------------------- |\n",
    "| tf.add(x, y, name=None)     | 求和                                                         |\n",
    "| tf.subtract(x, y, name=None)     | 减法                                                         |\n",
    "| tf.multiply(x, y, name=None)     | 乘法                                                         |\n",
    "| tf.divide(x, y, name=None)     | 除法                                                         |\n",
    "| tf.mod(x, y, name=None)     | 取模                                                         |\n",
    "| tf.abs(x, name=None)        | 求绝对值                                                     |\n",
    "| tf.neg(x, name=None)        | 取负 (y = -x).                                               |\n",
    "| tf.sign(x, name=None)       | 返回符号 y = sign(x) = -1 if x < 0; 0 if x == 0; 1 if x > 0. |\n",
    "| tf.inv(x, name=None)        | 取反                                                         |\n",
    "| tf.square(x, name=None)     | 计算平方 (y = x * x = x^2).                                  |\n",
    "| tf.round(x, name=None)      | 舍入最接近的整数 # ‘a’ is [0.9, 2.5, 2.3, -4.4] tf.round(a) ==> [ 1.0, 3.0, 2.0, -4.0 ] |\n",
    "| tf.sqrt(x, name=None)       | 开根号 (y = \\sqrt{x} = x^{1/2}).                             |\n",
    "| tf.pow(x, y, name=None)     | 幂次方  # tensor ‘x’ is [[2, 2], [3, 3]] # tensor ‘y’ is [[8, 16], [2, 3]] tf.pow(x, y) ==> [[256, 65536], [9, 27]] |\n",
    "| tf.exp(x, name=None)        | 计算e的次方                                                  |\n",
    "| tf.log(x, name=None)        | 计算log，一个输入计算e的ln，两输入以第二输入为底             |\n",
    "| tf.maximum(x, y, name=None) | 返回最大值 (x > y ? x : y)                                   |\n",
    "| tf.minimum(x, y, name=None) | 返回最小值 (x < y ? x : y)                                   |\n",
    "| tf.cos(x, name=None)        | 三角函数cosine                                               |\n",
    "| tf.sin(x, name=None)        | 三角函数sine                                                 |\n",
    "| tf.tan(x, name=None)        | 三角函数tan                                                  |\n",
    "| tf.atan(x, name=None)       | 三角函数ctan                                                 |\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "4b4c99d3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(), dtype=float32, numpy=2.0>"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = tf.constant(2.)\n",
    "y = tf.constant(3.)\n",
    "\n",
    "# 算术操作符：+ - * / % \n",
    "tf.add(x, y, name=None)        # 加法(支持 broadcasting)\n",
    "tf.subtract(x, y, name=None)   # 减法\n",
    "tf.multiply(x, y, name=None)   # 乘法\n",
    "tf.divide(x, y, name=None)     # 浮点除法, 返回浮点数(python3 除法)\n",
    "tf.math.mod(x, y, name=None)   # 取余  tf.mod(x, y, name=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "759ea985",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(), dtype=float32, numpy=0.6931472>"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 幂指对数操作符：^ ^2 ^0.5 e^ ln \n",
    "tf.pow(x, y, name=None)        # 幂次方\n",
    "tf.square(x, name=None)        # 平方\n",
    "tf.sqrt(x, name=None)          # 开根号，必须传入浮点数或复数\n",
    "tf.exp(x, name=None)           # 计算 e 的次方\n",
    "tf.math.log(x, name=None)      # 以 e 为底，必须传入浮点数或复数 tf.log(x, name=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "0d312d7a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(), dtype=float32, numpy=2.0>"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 取符号、负、倒数、绝对值、近似、两数中较大/小的\n",
    "tf.negative(x, name=None)      # 取负(y = -x).\n",
    "tf.sign(x, name=None)          # 返回 x 的符号\n",
    "tf.math.reciprocal(x, name=None)    # 取倒数\n",
    "tf.abs(x, name=None)           # 求绝对值\n",
    "tf.round(x, name=None)         # 四舍五入\n",
    "tf.math.ceil(x, name=None)     # 向上取整\n",
    "tf.floor(x, name=None)         # 向下取整\n",
    "tf.math.rint(x, name=None)     # 取最接近的整数 \n",
    "tf.maximum(x, y, name=None)    # 返回两tensor中的最大值 (x > y ? x : y)\n",
    "tf.minimum(x, y, name=None)    # 返回两tensor中的最小值 (x < y ? x : y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "6126f842",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(), dtype=float32, numpy=1.1071488>"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 三角函数和反三角函数\n",
    "tf.cos(x, name=None)    \n",
    "tf.sin(x, name=None)    \n",
    "tf.tan(x, name=None)    \n",
    "tf.acos(x, name=None)\n",
    "tf.asin(x, name=None)\n",
    "tf.atan(x, name=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "9ceecb66",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(), dtype=float32, numpy=1.0>"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 其它\n",
    "tf.divide(x, y, name=None)  # 计算 Python 风格的x除法y\n",
    "tf.math.divide_no_nan(x, y) # 计算一个安全除法，如果y（分母）为零，则返回 0\n",
    "tf.truediv(x, y, name=None) # Python 3 除法运算符语义\n",
    "tf.math.floordiv(x, y, name=None)  # 按元素除，向最负整数舍入\n",
    "tf.realdiv(x, y, name=None)   \n",
    "# tf.truncatediv(x, y, name=None)  # 某些函数被弃置了\n",
    "# tf.floor_div(x, y, name=None)\n",
    "tf.truncatemod(x, y, name=None)\n",
    "tf.math.floormod(x, y, name=None)  # 返回除法的元素余数\n",
    "# tf.cross(x, y, name=None)\n",
    "# tf.add_n(inputs, name=None)  # inputs: A list of Tensor objects, each with same shape and type\n",
    "tf.math.squared_difference(x, y, name=None) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5951da7",
   "metadata": {},
   "source": [
    "### 矩阵、张量运算\n",
    "\n",
    "https://www.tensorflow.org/api_docs/python/tf/linalg\n",
    "\n",
    "如果 AttributeError ，尝试  tf.linalg.xxx，或者 搜索 tf.linalg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "cdce9e27",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(2, 2), dtype=float32, numpy=\n",
       "array([[0.2       , 0.33333334],\n",
       "       [0.42857143, 0.5       ]], dtype=float32)>"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = tf.constant([[1., 2.], [3., 4.]])\n",
    "B = tf.constant([[5., 6.], [7., 8.]])\n",
    "\n",
    "tf.add(A, B)    # 计算矩阵A和B的和 C = A+B\n",
    "tf.subtract(A, B)    # 计算矩阵A和B的差 C = A-B\n",
    "tf.multiply(A, B)    # 计算矩阵A和B的乘积 C = A*B 哈达马积 Hadamard product\n",
    "tf.divide(A, B)    # 计算矩阵对应元素除法 C = A/B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "6b09a9bc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(), dtype=float32, numpy=10.0>"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "k = tf.constant(3.)\n",
    "tf.multiply(k, A)  # 数乘\n",
    "tf.matmul(A, B)    # 矩阵乘积，点积\n",
    "tf.reduce_sum(A, axis=None, keepdims=False, name=None) # 计算张量维度上的元素总和"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "d0274d1d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(), dtype=float32, numpy=5.4772253>"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tf.transpose(A, perm=None, name='transpose') # 转置，可以通过指定 perm=[1, 0] 来进行轴变换\n",
    "tf.linalg.trace(A, name=None)  # 计算张量的迹\n",
    "tf.linalg.det(A, name=None)    # 计算一个或多个方阵的行列式\n",
    "tf.linalg.inv(A, name=None)    # 计算一个或多个平方可逆矩阵或其伴随矩阵（共轭转置）的逆\n",
    "tf.linalg.svd(A, full_matrices=False, compute_uv=True, name=None) # 计算一个或多个矩阵的奇异值分解\n",
    "tf.linalg.qr(A, full_matrices=None, name=None)  # 计算一个或多个矩阵的 QR 分解\n",
    "tf.norm(A, ord='euclidean', axis=None, name=None) # 计算向量、矩阵和张量的范数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "416f7e9c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(3, 3), dtype=int32, numpy=\n",
       "array([[1, 0, 0],\n",
       "       [0, 2, 0],\n",
       "       [0, 0, 3]])>"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tf.eye(2, num_columns=None, batch_shape=None, dtype=tf.float32, name=None) # 构造一个单位矩阵，或一批矩阵\n",
    "tf.linalg.diag(diagonal=[1,2,3], name=None) # 构建一个对角矩阵。返回具有给定成批对角线值的成批对角线张量"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b25841b",
   "metadata": {},
   "source": [
    "## 导数"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f39fc696",
   "metadata": {},
   "source": [
    "### 一元显函数"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cab47bbb",
   "metadata": {},
   "source": [
    "**一阶导数**\n",
    "\n",
    "$\\large \\frac{d}{dx}(y = x^2)\\mid_{x=4}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "dd7140c6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tf.Tensor(16.0, shape=(), dtype=float32) tf.Tensor(8.0, shape=(), dtype=float32)\n"
     ]
    }
   ],
   "source": [
    "x = tf.Variable(initial_value=4.)\n",
    "with tf.GradientTape() as g:\n",
    "    y = tf.square(x)\n",
    "y_grad = g.gradient(y, x)\n",
    "print(y, y_grad)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f3d9528",
   "metadata": {},
   "source": [
    "$\\large \\frac{d}{dx}(y = x^2)\\mid_{x=3}$\n",
    "\n",
    "如果x是常数，需要加 g.watch(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "id": "7d201c1b",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tf.Tensor(6.0, shape=(), dtype=float32)\n"
     ]
    }
   ],
   "source": [
    "x = tf.constant(3.0)\n",
    "with tf.GradientTape() as g:\n",
    "    g.watch(x)\n",
    "    y = x * x\n",
    "dy_dx = g.gradient(y, x)\n",
    "print(dy_dx)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7474ad9d",
   "metadata": {},
   "source": [
    "$\\large \\frac{d}{dx}(z = y^2 = x^4)\\mid_{x=3}, \\frac{d}{dx}(y = x^2)\\mid_{x=3}, \\frac{d}{dy}(z = y^2)\\mid_{y=x^2=9}$\n",
    "\n",
    "$\\large 4x^3, 2x, 2y$\n",
    "\n",
    "设置persistent=True，可以多次求导"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "0229dc93",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tf.Tensor(108.0, shape=(), dtype=float32)\n",
      "tf.Tensor(6.0, shape=(), dtype=float32)\n",
      "tf.Tensor(18.0, shape=(), dtype=float32)\n"
     ]
    }
   ],
   "source": [
    "x = tf.constant(3.0)\n",
    "with tf.GradientTape(persistent=True) as g:\n",
    "    g.watch(x)\n",
    "    y = x * x\n",
    "    z = y * y\n",
    "dz_dx = g.gradient(z, x)  # (4*x^3 at x = 3)\n",
    "print(dz_dx)\n",
    "\n",
    "dy_dx = g.gradient(y, x)  # (2*x at x = 3)\n",
    "print(dy_dx)\n",
    "\n",
    "dz_dy = g.gradient(z, y)  # (2*y at y = x^2 = 9)\n",
    "print(dz_dy)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d358ffa",
   "metadata": {},
   "source": [
    "逐元素计算"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "id": "7476e761",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x233f55f2a00>"
      ]
     },
     "execution_count": 111,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEGCAYAAAB1iW6ZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAq5UlEQVR4nO3deXxU9b3/8dcnk4SENRDCDrIIsskSEJeruCGLC4jautZWa12qv9rb9lbutdp6rV3t5q3VUqu2dUFbl2JFwK2itci+BQQBWUJI2CGQdWa+vz/OgEOcwAAzOZPJ+/l45DEz53xn5jNnknfOfOec79ecc4iISOOX4XcBIiKSGAp0EZE0oUAXEUkTCnQRkTShQBcRSROZfj1x+/btXc+ePf16ehGRRmnhwoU7nHMFsdb5Fug9e/ZkwYIFfj29iEijZGYb61unLhcRkTShQBcRSRMKdBGRNOFbH3ostbW1FBcXU1VV5XcpDS4nJ4du3bqRlZXldyki0kilVKAXFxfTqlUrevbsiZn5XU6Dcc6xc+dOiouL6dWrl9/liEgjddQuFzN70sy2mdmKetabmT1iZmvNbJmZFR5vMVVVVeTn5zepMAcwM/Lz85vkJxMRSZx4+tCfBsYfYf0EoG/k51bgsRMpqKmF+UFN9XWLSOIctcvFOTfHzHoeockk4M/OG4d3rpnlmVln59zWRBUpIunJOUd1MEx1bZiqYIiaYJhg2BEKh6kNOUJhRzDsCIbCkUtHMByOXHrXQ2FH2Dmcg7DzHtM5cLjIbbz13hMeahN24Pis/cE2B68DhMOfDS/uDqs76jqHD0F++LrYK0b2bMfofjHPDTohiehD7wpsjrpdHFn2uUA3s1vx9uLp0aNHAp5aRPzinGNfZZDt+6vYtq+anQdq2FdVy77KIOVVtXWuB6moCVFVG/3jhXhTmpLh4Afx28/tk7KBHquvIOZb5JybCkwFGDlyZBN6G0UaH+ccJXur2LjzAJt2VrBxVwWbdlVQsqeSbfuq2b6/mppgOOZ9MzOM1rlZtM7JpHVuFq1yMmnXIpucrAA5mRneZdbBy8Ch29mBDDIDRmZGBpkZRiDDyApkEMiwz5YHjMyMz64HMowMMzIMDMPMC84Ms88uAaKuH1xnddocvE50m6jXFd01evjyw1+/X12oiQj0YqB71O1uQEkCHrfB3XfffbRv3567774bgHvvvZeOHTvyjW98w+fKRJIrGAqzams5y7fsZdXWfXxcuo+Pt5ZTXh081CYzw+jWNpeubXMZ1asdBa2a0aFVMwoiP+1bNqNNbhatc7LIycrQ90I+SESgTwfuMrNpwOnA3kT0nz/wWhErS/adcHHRBnZpzfcvG1Tv+q9+9atcccUV3H333YTDYaZNm8a8efMSWoNIKgiGwizZvIePPt3FvE93sXDjbvZHwrtls0z6d2rFpOFdOKVTa3q3b0GPds3p3CaHzIDORUxlRw10M3seOA9ob2bFwPeBLADn3OPADOBiYC1QAdyUrGKTrWfPnuTn57N48WLKysoYPnw4+fn5fpclkhCVNSHe+Xgbb60q452Pt7G3shaAvh1aMmlYF0b1akdhj7Z0a5urvetGKp6jXK49ynoH3JmwiiKOtCedTLfccgtPP/00paWl3Hzzzb7UIJIozjnmfbqLlxYVM2N5Kfurg+Q1z+LCAR0YM6AjZ/TOp12LbL/LlARJqTNFU8HkyZO5//77qa2t5bnnnvO7HJHjUlUbYvqSEp74YD1ryvbTIjvAxad2ZnJhV0b1bKeukzSlQK8jOzub888/n7y8PAKBgN/liByTqtoQf/n3Rn4/Zx079tfQv1Mrfn7VEC4Z0pnm2fpzT3d6h+sIh8PMnTuXv/71r36XIhK3UNjx8qJifvXmGkr2VnH2ye2547w+nNWn6Q2l0ZQp0KOsXLmSSy+9lMmTJ9O3b1+/yxGJy6qt+5jy0jKWFu9lSLc2PPyFoZx1cnu/yxIfKNCjDBw4kPXr1/tdhkhcqmpDPPruWh775zra5Gbx66uHMWlYF+2RN2EKdJFGaOPOA3z92UUUlezjisKu3HfJQNrqaJUmT4Eu0sjMLirl239digFP3DiSMQM7+l2SpAgFukgj4Zzj1299wm/e/oRTu7bhd9cX0r1dc7/LkhSiQBdpBEJhx/deXcHz8zZx1Yhu/PDyweRk6bBaOZzOLjiKH/zgBzz88MNHbPP888/z0EMPfW55z5492bFjR7JKkyaiqjbEnc8u4vl5m/j6eX34+VVDFOYSkwI9AWbOnMn48Uea1Enk+NQEw9z2l4XMLCrlvksH8t3x/XUUi9RLgR7DQw89xCmnnMKYMWNYvXo1oVCIwsLPpkr95JNPGDFiBOD1ay5ZsoTCwkJ27tzJ2LFjGT58OLfddhsuMnL//PnzGTJkCFVVVRw4cIBBgwaxYkXMKVpFDgmFHf/54hLeW7Odn1xxKl89WxOIy5Glbh/6G1OgdHliH7PTqTDhJ0dssnDhQqZNm8bixYsJBoMUFhYyYsQI2rRpw5IlSxg2bBhPPfUUX/nKVwBYvHgxQ4cOxcx44IEHOPvss7n//vt5/fXXmTp1KgCnnXYaEydO5Hvf+x6VlZXccMMNDB48OLGvTdKKc477/r6C15dt5X8u7s81ozTDlxxd6ga6T95//30mT55M8+be0QMTJ04EvFEYn3rqKX75y1/ywgsvHBonfebMmUyYMAGAOXPm8PLLLwNwySWX0LZt20OPe//993PaaaeRk5PDI4880pAvSRqhX731Cc995PWZ3zq6j9/lSCORuoF+lD3pZIrVR3nllVfywAMPcMEFFzBixIhD46TPnj2bl1566Yj3Bdi1axf79++ntraWqqoqWrRokZzipdGbuaKUR97+hC+O7MZ/jTvF73KkEVEfeh2jR4/mlVdeobKykvLycl577TUAcnJyGDduHHfccQc33eTN4bF3716CweChcB89ejTPPvssAG+88Qa7d+8+9Li33norDz74INdffz333HNPA78qaSzWbivn2y8uYWj3PB68fLC+AJVjokCvo7CwkKuvvpphw4Zx5ZVXcs455xxad/3112NmjB07FoA333yTMWPGHFr//e9/nzlz5lBYWMjs2bPp0cPr9/zzn/9MZmYm1113HVOmTGH+/Pm88847DfvCJOWVV9Vy618Wkpsd4PEbCmmWqUMT5djYwSMxGtrIkSPdggULDlu2atUqBgwY4Es98Xj44YfZu3cvDz74IOD1q99yyy2cccYZCXn8VH/9kjzOOe56bjEzi0p59pbTOaO3pj6U2MxsoXNuZKx1qduHnmImT57MunXrDtuzfuKJJ3ysSNLJ9KUlvL58K98df4rCXI6bAj1Or7zyit8lSJoq21fFfa+uYHiPPG7TES1yAlKuD92vLiC/NdXX3dQ555jy0jJqQmF+8YWhBDL0Jagcv5QK9JycHHbu3Nnkws05x86dO8nJyfG7FGlgLy7YzLurtzNlfH96F7T0uxxp5FKqy6Vbt24UFxezfft2v0tpcDk5OXTr1s3vMqQB7dxfzUOvr+KM3u248cyefpcjaSClAj0rK4tevTRehTQND89eTUVNiB9ePpgMdbVIAqRUl4tIU7GseA/T5m/mK2f15OQOrfwuR9KEAl2kgYXDjh9MLyK/RTO+Maav3+VIGlGgizSwVxZvYdGmPdwz/hRa52T5XY6kEQW6SAOqqg3xs1kfM6x7HlcW6ktwSSwFukgDembuRsr2VTNlQn99ESoJp0AXaSAHqoM89s91nH1ye53eL0mhQBdpIE9/uIGdB2r41th+fpciaSquQDez8Wa22szWmtmUGOvbmNlrZrbUzIrM7KbElyrSeO2trOX3763jwv4dKOzR9uh3EDkORw10MwsAjwITgIHAtWY2sE6zO4GVzrmhwHnAL8wsO8G1ijRaf3x/PfuqgvznRdo7l+SJZw99FLDWObfeOVcDTAMm1WnjgFbmTa/SEtgFBBNaqUgjVV5Vy1P/2sCEwZ0Y3LWN3+VIGosn0LsCm6NuF0eWRfstMAAoAZYDdzvnwnUfyMxuNbMFZragKY7XIk3T8/M2UV4d5I7zNDSuJFc8gR7r2Kq6wyGOA5YAXYBhwG/NrPXn7uTcVOfcSOfcyIKCgmMsVaTxqQmGefKDDZzZO58h3fL8LkfSXDyBXgx0j7rdDW9PPNpNwMvOsxb4FOifmBJFGq/pS0so3VfFref29rsUaQLiCfT5QF8z6xX5ovMaYHqdNpuACwHMrCNwCrA+kYWKNDbOOf4wZz2ndGzFef30iVSS76iB7pwLAncBs4BVwIvOuSIzu93Mbo80exA4y8yWA28D9zjndiSraJHG4J9rtrO6rJyvje6Nd7yASHLFNR66c24GMKPOssejrpcAYxNbmkjj9oc56+nUOoeJQ7v4XYo0ETpTVCQJ1pSV8+G6ndx41klkZ+rPTBqGftNEkuDZuRvJDmRw9cjuR28skiAKdJEEO1Ad5OVFW5hwaifyWzbzuxxpQhToIgk2fWkJ5dVBvnTGSX6XIk2MAl0kgZxzPDN3I/07tWLESRqESxqWAl0kgZZs3kNRyT6uP+MkHaooDU6BLpJAz8zdRIvsAJOH1x3uSCT5FOgiCbKvqpZ/LCvh8uFdadksrlM8RBJKgS6SIK8v20p1MMzVp+lQRfGHAl0kQf62sJi+HVpyqsY8F58o0EUSYP32/SzcuJurRnTTl6HiGwW6SAK8tKiYDENfhoqvFOgiJygUdry8aAvn9iugQ+scv8uRJkyBLnKCPly3g617q7hqhL4MFX8p0EVO0N8WFtMmN4sLB3TwuxRp4hToIidgf3WQWUWlXDa0MzlZAb/LkSZOgS5yAt5cWUpVbZjLh+nLUPGfAl3kBExfUkLXvFwKe2ggLvGfAl3kOO0+UMP7n+zg0qGdycjQsefiPwW6yHGasWIrwbDTnKGSMhToIsdp+pIS+hS0YGDn1n6XIgIo0EWOS+neKuZt2MXEoV11qr+kDAW6yHH4x7ISnIOJw9TdIqlDgS5yHKYvLeHUrm3o1b6F36WIHKJAFzlGm3dVsKx4L5cN7ex3KSKHUaCLHKNZRaUAjB+kQJfUokAXOUazV5bRv1MreuQ397sUkcMo0EWOwc791SzYsIuxgzr5XYrI5yjQRY7B26u2EXYwblBHv0sR+RwFusgxmFVUSte8XJ1MJCkprkA3s/FmttrM1prZlHranGdmS8ysyMzeS2yZIv7bXx3k/bU7GDeok04mkpSUebQGZhYAHgUuAoqB+WY23Tm3MqpNHvA7YLxzbpOZaaR/STtz1mynJhhWd4ukrHj20EcBa51z651zNcA0YFKdNtcBLzvnNgE457YltkwR/80qKqVdi2xG9mzndykiMcUT6F2BzVG3iyPLovUD2prZP81soZndGOuBzOxWM1tgZgu2b99+fBWL+KAmGOadj7cxZkAHAhoqV1JUPIEe67fX1bmdCYwALgHGAfeZWb/P3cm5qc65kc65kQUFBcdcrIhf5q7fSXlVkLEDdbiipK6j9qHj7ZFHT2feDSiJ0WaHc+4AcMDM5gBDgTUJqVLEZ7OKSmmeHeDsvu39LkWkXvHsoc8H+ppZLzPLBq4Bptdp83fgHDPLNLPmwOnAqsSWKuKPcNjx5soyzu1XoImgJaUddQ/dORc0s7uAWUAAeNI5V2Rmt0fWP+6cW2VmM4FlQBh4wjm3IpmFizSUJcV72FZezTidHSopLp4uF5xzM4AZdZY9Xuf2z4GfJ640kdQwu6iMzAzj/P46GldSm84UFTkC5xyzi0o5s08+bXKz/C5H5IgU6CJHsG77ftbvOKDBuKRRUKCLHMGsojIALhqgs0Ml9SnQRY5gVlEpw7rn0alNjt+liByVAl2kHiV7KllWvJexGrtFGgkFukg93lzpdbfocEVpLBToIvWYvbKUkzu0pE9BS79LEYmLAl0khj0VNcxdv4uxA9XdIo2HAl0khrdXbSMUdupukUZFgS4Sw+yVpXRqncOpXdv4XYpI3BToInVU1oR4b812xg7qSIbGPpdGRIEuUsf7n2ynqjassc+l0VGgi9Qxq6iMNrlZnN5bU81J46JAF4kSDIV5++MyLuzfgayA/jykcdFvrEiUeRt2saeiVmeHSqOkQBeJMruojGaZGYzupzlvpfFRoItEHBz7/Jy+BTTPjmvuF5GUokAXiVixZR8le6sYp+4WaaQU6CIRs4pKyTC4UGOfSyOlQBeJmL2ylFG92tGuRbbfpYgcFwW6CPDpjgOsKduvsVukUVOgiwCzi0oBuEijK0ojpkAXwes/H9y1Nd3aNve7FJHjpkCXJm/bvioWb96jsVuk0VOgS5P35qoynNNUc9L4KdClyZtVVMZJ+c3p11FTzUnjpkCXJm1fVS3/XreDcYM6Yaaxz6VxU6BLk/bux9uoDTnNHSppQYEuTdrsojLat2xGYY+2fpcicsIU6NJkVdWGeHf1Nk01J2lDgS5N1ntrtlNRE2LCYB3dIukhrkA3s/FmttrM1prZlCO0O83MQmZ2VeJKFEmOmStKaZObxRm98/0uRSQhjhroZhYAHgUmAAOBa81sYD3tfgrMSnSRIolWEwzz1qoyLhrYUVPNSdqI5zd5FLDWObfeOVcDTAMmxWj3/4CXgG0JrE8kKf61bgflVUF1t0haiSfQuwKbo24XR5YdYmZdgcnA40d6IDO71cwWmNmC7du3H2utIgkzc3kpLZtlcnbf9n6XIpIw8QR6rK//XZ3bvwbucc6FjvRAzrmpzrmRzrmRBQWas1H8EQyFmb2ylAv6d6BZZsDvckQSJp6JE4uB7lG3uwElddqMBKZFzrRrD1xsZkHn3KuJKFIkkeZ9uovdFbXqbpG0E0+gzwf6mlkvYAtwDXBddAPnXK+D183saeAfCnNJVW+sKCUnK4NzT9GnREkvRw1051zQzO7CO3olADzpnCsys9sj64/Yby6SSsJhx6yiUs7r14Hm2fHsz4g0HnH9RjvnZgAz6iyLGeTOua+ceFkiybFo0262lVcz4VR1t0j60QG40qS8saKU7EAGF/Tv4HcpIgmnQJcmwznHzBWlnN23Pa1ysvwuRyThFOjSZCzfspcteyoZr6NbJE0p0KXJeG1pCVkB09jnkrYU6NIkhMOOfyzbyui+BeQ1z/a7HJGkUKBLkzB/wy627q1i4rAufpcikjQKdGkSpi8tIScrgzED1N0i6UuBLmmvNhRmxvKtjBnQkRbNdDKRpC8FuqS9D9buYHdFLROHqrtF0psCXdLea0tLaJWTqbFbJO0p0CWtVdWGmF1UxvhBnTRUrqQ9BbqktbdXbWN/dVBHt0iToECXtPbSomI6tm7GWX00M5GkPwW6pK1t+6p4b812rijsRiAj1sRbIulFgS5p69UlWwiFHVcWdvO7FJEGoUCXtOSc46WFWxjeI4+TO7T0uxyRBqFAl7S0Yss+VpeVc9UI7Z1L06FAl7T0t4Wbyc7M4NIhOrpFmg4FuqSd6mCIvy8tYdygTrTJ1UQW0nQo0CXtvLVyG3sqatXdIk2OAl3SzrMfbaRrXi5nn6xjz6VpUaBLWlm3fT8frtvJdaf30LHn0uRoLFFJK8/O3URWwPjiyO7HfufaSti9Ear2QHZLyOsOOW0SXqNIsijQJW1U1oT428LNjBvUiYJWzeK7U/V+WPE3WPoCFM+HcG3USoOOg2HQ5TD8S9BKk2NIalOgS9p4bWkJ+6qCfOmMk47eOBSEBU/Cez+Fih1QMADO/Dp0GgK5eVBzALavhnXvwjsPwpyfw+m3wznfhpzWSX8tIsdDgS5p45mPNtKvY0tG9Wp35IY71sKrt3t75D3PgQu+B91PB4vR537ud732c34O//oNrHgZLv8d9DonOS9C5AToS1FJC8uK97CseC/Xn34SFiuYD/p4Bkw9D3Z8Alf+Eb78GvQ4I3aYH9T+ZLji93DzLAhkwp8neuHuXMJfh8iJUKBLWvjjB5/SIjvA5MKu9Tea+xhMuxby+8Ad/4JTrzpykNfV43S47X0YMBHevB9euxvC4RMvXiRBFOjS6BXvruAfy7ZyzagetM6p58zQd38EM6fAgMvg5pnQ5jhPOmrWEr7wtNeXvuhP8OodXn+8SApQH7o0ek9+sAGAm8/uFbvBnIe9Lz+H3wCXPQIZJzgVnRlceD9k5sK7P4RQNVzxBwhomAHxV1x76GY23sxWm9laM5sSY/31ZrYs8vOhmQ1NfKkin7e3opZp8zcxcWgXuublfr7B3Me9o1SGXA2X/d+Jh3m0c/8Lxv4Qil6BV25T94v47qh76GYWAB4FLgKKgflmNt05tzKq2afAuc653WY2AZgKnJ6MgkWiPfPRRipqQnztnN6fX7n4GZh5D/S/FCb9DjKS0MN41v+DcBDe+gG06gzjHkr8c4jEKZ4ul1HAWufcegAzmwZMAg4FunPuw6j2cwGNiiRJVx0M8fSHGzinb3sGdqlzbPinc7wvLXufD1c96R2dkiz/8U3YVwL//i207gJn3pm85xI5gnh2WboCm6NuF0eW1eerwBuxVpjZrWa2wMwWbN++Pf4qRWJ4edEWtpdXc9voPoev2PUpvPhlaNcHvvhnyIzzrNHjZQbjf+J94Trrf7wuGBEfxBPosY7rinkArpmdjxfo98Ra75yb6pwb6ZwbWVBQEH+VInVUB0P89p21DO2ex3+cnB+1ohymXQcuDNc+33BndWYEvC9Gu58Or9wBW5c2zPOKRIkn0IuB6JGOugEldRuZ2RDgCWCSc25nYsoTie2F+ZvZsqeS74zt99mJROEwvHybd8r+F//kHW/ekLJy4epnoHk7eP462K9PodKw4gn0+UBfM+tlZtnANcD06AZm1gN4GfiSc25N4ssU+UxlTYj/e2cto3q1O3zM83/+CFa/DuN/DL3P86e4lh3gmme98WFevBGCNf7UIU3SUQPdORcE7gJmAauAF51zRWZ2u5ndHml2P5AP/M7MlpjZgqRVLE3eM3M3sr28mu+MPeWzvfMVL3njrQz/Eoy61d8CuwyHSY/Cpg/hje/6W4s0KXF99e+cmwHMqLPs8ajrtwC3JLY0kc/bXx3ksffWMbpfwWeDcJUsgVfvhO5nwCW/OLbT+ZPl1KugdDn869fQaTCcpj8PST6d+i+Nyh/mrGfXgRq+fVE/b8H+bTDtemieD1f/JflHtByLC++HvuPgjXtgwwd+VyNNgAJdGo3Nuyp4/L11XDKkM0O750GwGl64ASp2wrXPef3XqSQjAFf+Adr19vrTd2/0uyJJcwp0aTR+NGMVGWbce/EAb+ja178Fmz+CyY9B5xQdbSKnDVzzvDeA17TrvBmSRJJEgS6Nwr/W7uCNFaXceX4fuuTlwke/907tH/1dGDTZ7/KOrP3J8IUnYdtKb3RGjfkiSaJAl5RXGwrz/elF9GjXnFvO6Q2fvAWz/tsbo+W8//a7vPicPAYu+l9YNd07GkckCTR8rqS8p/71KWu37eeJG0eSs/sT+NtN0GEQTP59cgbcSpYz74KyIu94+Y4DvaECRBKoEf01SFO0dls5D89ew5gBHbmwRwY890XvjMzrpnmTTTQmZnDpr6HrSO+M1rIivyuSNKNAl5QVDIX59otLaZEd4EcT+2Iv3AD7y7wvGY93xiG/ZeV4wwM0awXPXwMHdvhdkaQRBbqkrMffW8fS4r38cNJgOrz7Xdg8Fy5/DLqN8Lu0E9O6M1zznHcM/XNXQ02F3xVJmlCgS0oqKtnLb97+hMuGduGSPc/Asmlw/r0w+Aq/S0uMbiPgyj9CySJ46aual1QSQoEuKWdfVS13PbeYvObZ/KTHfHj3IRh6LYz+L79LS6wBl8KEn8HqGd6YLy7mqNQicdNRLpJSwmHHt15YwuZdFcy8aAct3vyud/r8xP9LjTFaEm3U12BvsTfmS/N8uOBevyuSRkyBLinlt++u5a1V25h61h5Ofv9b0OMM+MLTEMjyu7TkufD73vAFc37mfWl6zrf9rkgaKQW6pIx3Pi7jV2+tYUrfLVy09HtQ0B+unQbZzf0uLbkyMuCy30CwCt7+X8jMhTO/7ndV0ggp0CUlLNy4mzufXcyX8tdwW8mPsYJ+cON0yM3zu7SGkRGAyx/3BhybFTn7VaEux0hfiorvVm3dx01PzWNy86U8UPkjrEN/L8ybt/O7tIYVyPSOfBlwmRfq7/5IX5TKMVGgi6827jzAl/44j+sDb/NQzU+wjoPhxr83vTA/KDMbrnoaht0A7/3UG0tdg3lJnNTlIr5ZU1bOjU/8m68Hn+VmXoW+Y+GqpxrfKf2JFsj0jurJaQNzH4XyEm/cmuwWflcmKU576OKLhRt3ceNj7/Bg8FdemI/4indKf1MP84MyMmDcQzD2Ifj4dfjjONizye+qJMUp0KXBvbWyjPueeJlp9j+MYS6MecAbtCqgD4yHMYOz7oLrXoQ9G2Hq+fDpHL+rkhSmQJcGEwo7Hp75Ma8982teCtxLj5wq7EuvwtnfTM+ThhKl70Vwy9uQ2xb+NBHeegBCtX5XJSlIgS4NYsf+au6cOotB/7qL32T/jmbdhpBxx/vQ+1y/S2scCvrBbe/B8Bvgg1/Ck+Ngx1q/q5IUo8+4klTOOV5bWsJHf/89Pwo/SZusKrjgATLOvEtdLMcquwVM+i2cfCG8djc8dqZ3VunZ/wmZzfyuTlKA/qIkaUr2VPLECy8xYctveChjDZUdhhL4wlTo0N/v0hq3QZOhx1kw63/gnz+G5X+FcT/2umbUddWkKdAl4fZW1DJt5rt0XPII37MPqM5pS3jcI+QOv8E7I1JOXKuOcNUfYdh1MOM78NwXvJAf8wPocbrf1YlPzPl0JtrIkSPdggULfHluSY69FbXMeOdt2ix4hHHuQ0IZ2VQX3kKri6ZATmu/y0tfwRpY/Gd472fejE69z/eOjulzofbY05CZLXTOjYy5ToEuJ2pD2R7mzfwLJ61/jtNtJVWWw/4hN9H+om9DywK/y2s6ag7AvKkw93HYXwodBnrD8w6+0jtJSdKCAl0Sbl9lNfPem0lw2V8ZceA9CmwfO7M6ERx+Mx3P+1rTPXU/FQRrYMXf4N+PQtkKyMyBARNhyBeh12h9gdrIHSnQ1YcucSvZvoPV/34dt2YWA8r/zRjbRTXZbC44h8yzbyJ/yMXqI08Fmdle3/rQa6FkMSx+xgv45S9Cs9beEAv9xkHPc7z5TSVtaA9dYnLOsXnzRjYveYfaDf+mw57F9A2tJ8tCVJDDprzTyRkyiZPOugrTx/nUF6yG9e/BqunelHcVO73l+X2h1znQ82zoUghte6rfPcWpy0Xq5Zxj9+5dlG1cxZ4NSwltLaL53jV0qV5PJ7w/+mqy2JTTn4pOp9H+1IvoMuQCLCvH58rluIVDULocNrwPn74PGz+EmnJvXbM20OlU6DwEOg6Cdn0gvw+0KFDQp4gTDnQzGw/8BggATzjnflJnvUXWXwxUAF9xzi060mMq0JMvHAqzd88O9pRtpnznFqp2lRDcVwr7y8g6sJXWlcUUBEtpZ+WH7lPjMtmS1YM9rfriOgyiw6Bz6TrgDAV4OgsFoWw5bF0KW5d5l2VFEKz8rE2z1tCuN7TpBq06QavO3k/ryGXzfMjJ87p7JKlOqA/dzALAo8BFQDEw38ymO+dWRjWbAPSN/JwOPBa5FCAcChEKBQkFawkGawkFg4SCNYTDIUKR2+FgLaFQkHColnCwltqaSoLVFYRqqgjXVBCqriRcW4mrrcJFLglWYbUVBGrKyawtJztYTrPQAZqH99PC7aeFq6SthWlbp55ql8XOjHx2NevCuryBrG3bk2YFfSjoPYROPQfRK0t/lE1KIBO6DPd+DgoFvQHBdq2Hnetg17rI5XrY8AFU7Yn9WFktvFmmcvK8y9y23j+D7OaQlQtZBy9bRC4PLsuBjCxv7thA1mfXM7K8+g7dzvxseUYALCPyo08PEN+XoqOAtc659QBmNg2YBEQH+iTgz87b3Z9rZnlm1tk5tzXRBS9792+0ef/7ABgOi/qEYTjARS7Be4s/a3PodlSbQ49Tz+3o+3y27vOPG+sxAoQJECJgjgwg0dMc17oAVZbNAVpQmdGSqkALDmQXsCerN+FmrQk3a43ltiMrrzM5bbvQuqArbTt0p0XrdnQxo0uC65E0Esj0ulry+3hnoNZVU+EdGrlvK5RvhcrdULnHC/rKPd7tqj2w61Oo3ge1FVBb6V0mi2UAFhXydX8s8hNjHRb1TyHqn4PVuXLYP466y6LvZ/W0iSwrvNE7VyDB4gn0rsDmqNvFfH7vO1abrsBhgW5mtwK3AvTo0eNYawUgu2UeO5v3ORSr3oayOrfhUBRH1kdv2JhtP7c8+s3J+Kxd9GPY4Y9b9z4uI+DtUWRkensTGVlYZNmhy4C3PiOQCYFMMiLtA9k5ZGbnEmiWS1azFmTn5JLdrDnZuS3IzmlOTm4LsjKzyAJaHdeWFDkB2c29Lph2vY/tfuGwNxn2wXCPvgzVQLjW+3QQrvVGlAwHI5cxbruwN0WfC9f5qbuszm3quQ/UmfKv7rKodXWXHev9WnY4tu0Wp3gCPdZnmbod7/G0wTk3FZgKXh96HM/9Of1PGwOnjTmeu4qI3zIyvH8G2c2BfL+rSTvxDJ9bDHSPut0NKDmONiIikkTxBPp8oK+Z9TKzbOAaYHqdNtOBG81zBrA3Gf3nIiJSv6N2uTjngmZ2FzAL77DFJ51zRWZ2e2T948AMvEMW1+IdtnhT8koWEZFY4jr13zk3Ay+0o5c9HnXdAXcmtjQRETkWmoJORCRNKNBFRNKEAl1EJE0o0EVE0oRvoy2a2XZg43HevT2wI4HlJEqq1gWpW5vqOjaq69ikY10nOediTgXmW6CfCDNbUN9oY35K1bogdWtTXcdGdR2bplaXulxERNKEAl1EJE001kCf6ncB9UjVuiB1a1Ndx0Z1HZsmVVej7EMXEZHPa6x76CIiUocCXUQkTaRsoJvZF8ysyMzCZjayzrr/NrO1ZrbazMbVc/92ZvammX0Suaw7tWYianzBzJZEfjaY2ZJ62m0ws+WRdkmfGdvMfmBmW6Jqu7ieduMj23CtmU1pgLp+bmYfm9kyM3vFzPLqadcg2+torz8yHPQjkfXLzKwwWbVEPWd3M3vXzFZFfv/vjtHmPDPbG/X+3p/suqKe+4jvjU/b7JSobbHEzPaZ2TfrtGmQbWZmT5rZNjNbEbUsrixKyN+jcy4lf4ABwCnAP4GRUcsHAkuBZkAvYB0QiHH/nwFTItenAD9Ncr2/AO6vZ90GoH0DbrsfAN85SptAZNv1BrIj23RgkusaC2RGrv+0vvekIbZXPK8fb0joN/Bm5DoD+KgB3rvOQGHkeitgTYy6zgP+0VC/T8fy3vixzWK8r6V4J980+DYDRgOFwIqoZUfNokT9PabsHrpzbpVzbnWMVZOAac65aufcp3hjsI+qp92fItf/BFyelELx9kqALwLPJ+s5kuDQ5N/OuRrg4OTfSeOcm+2cC0ZuzsWb2cov8bz+Q5OfO+fmAnlm1jmZRTnntjrnFkWulwOr8ObnbSwafJvVcSGwzjl3vGehnxDn3BxgV53F8WRRQv4eUzbQj6C+Canr6ugisyZFLpMzK6vnHKDMOfdJPesdMNvMFkYmym4Id0U+8j5Zz0e8eLdjstyMtycXS0Nsr3hev6/byMx6AsOBj2KsPtPMlprZG2Y2qKFq4ujvjd+/V9dQ/46VX9ssnixKyHaLa4KLZDGzt4BOMVbd65z7e313i7EsacdexlnjtRx57/w/nHMlZtYBeNPMPo78J09KXcBjwIN42+VBvO6gm+s+RIz7nvB2jGd7mdm9QBB4tp6HSfj2ilVqjGXHNfl5MphZS+Al4JvOuX11Vi/C61LYH/l+5FWgb0PUxdHfGz+3WTYwEfjvGKv93GbxSMh28zXQnXNjjuNu8U5IXWZmnZ1zWyMf+bYlo0YzywSuAEYc4TFKIpfbzOwVvI9XJxRQ8W47M/sD8I8Yq5IysXcc2+vLwKXAhS7SeRjjMRK+vWJI2cnPzSwLL8yfdc69XHd9dMA752aY2e/MrL1zLumDUMXx3vg5YfwEYJFzrqzuCj+3GfFlUUK2W2PscpkOXGNmzcysF95/2Xn1tPty5PqXgfr2+E/UGOBj51xxrJVm1sLMWh28jvfF4IpYbROlTp/l5HqeL57JvxNd13jgHmCic66injYNtb1ScvLzyPcxfwRWOed+WU+bTpF2mNkovL/jncmsK/Jc8bw3fk4YX+8nZb+2WUQ8WZSYv8dkf+t7vD94QVQMVANlwKyodffifSO8GpgQtfwJIkfEAPnA28Ankct2SarzaeD2Osu6ADMi13vjfWO9FCjC63pI9rb7C7AcWBb5pehct67I7YvxjqJY10B1rcXrJ1wS+Xncz+0V6/UDtx98P/E+Bj8aWb+cqKOtkljT2XgftZdFbaeL69R1V2TbLMX7cvmsZNd1pPfG720Wed7meAHdJmpZg28zvH8oW4HaSH59tb4sSsbfo079FxFJE42xy0VERGJQoIuIpAkFuohImlCgi4ikCQW6iEiaUKCLiKQJBbqISJpQoItEmNlpkQHNciJnRRaZ2WC/6xKJl04sEoliZj8EcoBcoNg592OfSxKJmwJdJEpkHI35QBXe6eEhn0sSiZu6XEQO1w5oiTdbUI7PtYgcE+2hi0Qxs+l4s8X0whvU7C6fSxKJm6/joYukEjO7EQg6554zswDwoZld4Jx7x+/aROKhPXQRkTShPnQRkTShQBcRSRMKdBGRNKFAFxFJEwp0EZE0oUAXEUkTCnQRkTTx/wGw/D/M69kLuAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = tf.linspace(-10.0, 10.0, 200+1)\n",
    "\n",
    "with tf.GradientTape() as t:\n",
    "    t.watch(x)\n",
    "    y = tf.nn.sigmoid(x)\n",
    "\n",
    "dy_dx = t.gradient(y, x)\n",
    "\n",
    "plt.plot(x, y, label='y')\n",
    "plt.plot(x, dy_dx, label='dy/dx')\n",
    "plt.xlabel('x')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "647d4a89",
   "metadata": {},
   "source": [
    "**二阶导数**\n",
    "\n",
    "$\\large \\frac{d}{dx}(y = x^2)\\mid_{x=5}, \\frac{d^2}{dx^2}(y = x^2)\\mid_{x=5}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "668b90b7",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tf.Tensor(10.0, shape=(), dtype=float32)\n",
      "tf.Tensor(2.0, shape=(), dtype=float32)\n"
     ]
    }
   ],
   "source": [
    "x = tf.constant(5.0)\n",
    "with tf.GradientTape() as g:\n",
    "    g.watch(x)\n",
    "    with tf.GradientTape() as gg:\n",
    "        gg.watch(x)\n",
    "        y = x * x\n",
    "    dy_dx = gg.gradient(y, x)   # dy_dx = 2 * x\n",
    "d2y_dx2 = g.gradient(dy_dx, x)  # d2y_dx2 = 2\n",
    "print(dy_dx)\n",
    "\n",
    "print(d2y_dx2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f95c7a7b",
   "metadata": {},
   "source": [
    "**n阶导数**\n",
    "\n",
    "$\\large \\frac{d^n}{dx^n}(y = f(x))\\mid_{x=x_0}$\n",
    "\n",
    "需要重写，我不精通数据结构与算法，写不出像 mathematica 那样的用法。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "id": "8980657d",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def diff(func:str, x, n:int, ifexec=False):\n",
    "    \"\"\"\n",
    "    func: a string, such as func=\"y = x * x\";\n",
    "    x: tf.FloatTensor, such as x=tf.constant(1.0);\n",
    "    n: a int number, Order of derivative, such as n=1;\n",
    "    ifexec: \n",
    "    return dny_dxn: the value of order of derivative.\n",
    "    example: \n",
    "        x = tf.constant(1.0)\n",
    "        n = 3\n",
    "        func = \"y = x * x * x\"\n",
    "        diff(func, x, n)\n",
    "    \"\"\"\n",
    "    if n<1:\n",
    "        global y\n",
    "        exec(func, globals())\n",
    "    else:\n",
    "        with tf.GradientTape() as g:\n",
    "            g.watch(x)\n",
    "            diff(func, x, n-1,ifexec)\n",
    "            try:\n",
    "                y = g.gradient(y, x)\n",
    "            except TypeError:\n",
    "                pass\n",
    "    # print(n, y)\n",
    "    if y== None: y = tf.constant(0.)\n",
    "    if ifexec:\n",
    "        # exec executes a Python statement stored in a string or file\n",
    "        exec(\"d{0}y_dx{0} = y\".format(n), globals())\n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "63b99f52",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(), dtype=float32, numpy=6.0>"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# diff(func=\"y = x * x * x\", x=tf.constant(1.0), n=3)  \n",
    "# 这是一个错误用法，但是我还不知道为什么错\n",
    "\n",
    "x = tf.constant(1.0)\n",
    "diff(func=\"y = x * x * x\", x=x, n=3)    # 这是正确用法"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4714acf",
   "metadata": {},
   "source": [
    "$\\large \\frac{d^5}{dx^5}(y = 2sin(x))\\mid_{x=0}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "2d1438f6",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(), dtype=float32, numpy=2.0>"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = tf.constant(0.0)\n",
    "n = 5\n",
    "func = \"y = 2*tf.sin(x)\"\n",
    "diff(func, x, n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "id": "112edec8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "d0y_dx0:  tf.Tensor(1.0, shape=(), dtype=float32)\n",
      "d1y_dx1:  tf.Tensor(4.0, shape=(), dtype=float32)\n",
      "d2y_dx2:  tf.Tensor(12.0, shape=(), dtype=float32)\n",
      "d3y_dx3:  tf.Tensor(24.0, shape=(), dtype=float32)\n",
      "d4y_dx4:  tf.Tensor(24.0, shape=(), dtype=float32)\n",
      "d5y_dx5:  tf.Tensor(0.0, shape=(), dtype=float32)\n",
      "d6y_dx6:  tf.Tensor(0.0, shape=(), dtype=float32)\n"
     ]
    }
   ],
   "source": [
    "x = tf.constant(1.)\n",
    "n = 7\n",
    "func = \"y = x * x * x * x\"\n",
    "diff(func, x, n, ifexec=True)\n",
    "\n",
    "for i in range(n):\n",
    "    dny_dxn = \"print('d{0}y_dx{0}: ', d{0}y_dx{0})\".format(i)\n",
    "    exec(dny_dxn, globals())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "565bc494",
   "metadata": {},
   "source": [
    "**误差**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "id": "05a93a32",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(), dtype=float32, numpy=72.0>"
      ]
     },
     "execution_count": 132,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = tf.constant(1.)\n",
    "# x = tf.constant(100000.)\n",
    "n = 2\n",
    "func = \"y = x * x * x * x * x * x * x * x * x\"  # 误差更小\n",
    "# func = \"y = tf.pow(x,9)\"          # x增大时，tf.pow 会有误差,等价于 \"y = x**9\"\n",
    "diff(func, x, n)    # y in [-3.4E+38,3.4E+38]\n",
    "\n",
    "# for i in range(n):\n",
    "#     dny_dxn = \"print('d{0}y_dx{0}: ', d{0}y_dx{0})\".format(i)\n",
    "#     exec(dny_dxn, globals())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "id": "63c1e31a",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(<tf.Tensor: shape=(), dtype=float32, numpy=100000010.0>,\n",
       " <tf.Tensor: shape=(), dtype=int32, numpy=100000000>)"
      ]
     },
     "execution_count": 125,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tf.pow(100.,4), tf.pow(100,4)  # 浮点数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "id": "8202482b",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(<tf.Tensor: shape=(), dtype=float32, numpy=9.999981e+19>,\n",
       " <tf.Tensor: shape=(), dtype=int32, numpy=1661992960>)"
      ]
     },
     "execution_count": 124,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tf.pow(100.,10), tf.pow(100,10) # 溢出"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0537a4ea",
   "metadata": {},
   "source": [
    "**基本求导公式表**的实现："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2eaa1026",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "c504c335",
   "metadata": {},
   "source": [
    "## 偏导数\n",
    "\n",
    "https://www.tensorflow.org/guide/autodiff\n",
    "\n",
    "https://www.tensorflow.org/guide/advanced_autodiff"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e4feeff",
   "metadata": {},
   "source": [
    "### 多元显函数"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "457ae661",
   "metadata": {},
   "source": [
    "**一阶偏导**\n",
    "\n",
    "2元函数\n",
    "\n",
    "标量\n",
    "\n",
    "$\\large z = (x+2y)^2 = x^2+4y^2+4xy$\n",
    "\n",
    "$\\large \\frac{\\partial}{\\partial x}(z = (x+2y)^2)\\mid_{(x,y)=(1,0)}, \\frac{\\partial}{\\partial y}(z = (x+2y)^2)\\mid_{(x,y)=(1,0)}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "id": "40b1a1b5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tf.Tensor(2.0, shape=(), dtype=float32) tf.Tensor(4.0, shape=(), dtype=float32)\n"
     ]
    }
   ],
   "source": [
    "x = tf.constant(1.0)\n",
    "y = tf.constant(0.0)\n",
    "with tf.GradientTape(persistent=True) as g:\n",
    "    g.watch(x)\n",
    "    g.watch(y)\n",
    "    z = (x+2*y)*(x+2*y)\n",
    "    \n",
    "dz_dx = g.gradient(z, x)\n",
    "dz_dy = g.gradient(z, y)\n",
    "print(dz_dx, dz_dy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "2f4b3a0a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tf.Tensor(12.0, shape=(), dtype=float32) tf.Tensor(24.0, shape=(), dtype=float32) tf.Tensor(36.0, shape=(), dtype=float32)\n"
     ]
    }
   ],
   "source": [
    "x = tf.constant(1.0)\n",
    "y = tf.constant(1.0)\n",
    "m = tf.constant(1.0)\n",
    "with tf.GradientTape(persistent=True) as g:\n",
    "    g.watch(x)\n",
    "    g.watch(y)\n",
    "    g.watch(m)\n",
    "    temp = 1*x\n",
    "    temp += 2*y\n",
    "    temp += 3*m\n",
    "    z = tf.pow(temp,2)\n",
    "    \n",
    "dz_dx = g.gradient(z, x)\n",
    "dz_dy = g.gradient(z, y)\n",
    "dz_dm = g.gradient(z, m)\n",
    "print(dz_dx, dz_dy, dz_dm)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f41f7ffe",
   "metadata": {},
   "source": [
    "向量，多元函数\n",
    "\n",
    "$\\large z = (x_1+2x_2+\\dots+nx_n)^2 = (\\vec{k} \\vec{x})^2$\n",
    "\n",
    "$\\large \\frac{\\partial}{\\partial x_i}(z)\\mid_{\\vec{x}=\\vec{1}}$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26d7d5c3",
   "metadata": {},
   "source": [
    "第一种方法：通过 exec 自动命名变量实现"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "id": "c75c7149",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dz_dx1:  tf.Tensor(110.0, shape=(), dtype=float32)\n",
      "dz_dx2:  tf.Tensor(220.0, shape=(), dtype=float32)\n",
      "dz_dx3:  tf.Tensor(330.0, shape=(), dtype=float32)\n",
      "dz_dx4:  tf.Tensor(440.0, shape=(), dtype=float32)\n",
      "dz_dx5:  tf.Tensor(550.0, shape=(), dtype=float32)\n",
      "dz_dx6:  tf.Tensor(660.0, shape=(), dtype=float32)\n",
      "dz_dx7:  tf.Tensor(770.0, shape=(), dtype=float32)\n",
      "dz_dx8:  tf.Tensor(880.0, shape=(), dtype=float32)\n",
      "dz_dx9:  tf.Tensor(990.0, shape=(), dtype=float32)\n",
      "dz_dx10:  tf.Tensor(1100.0, shape=(), dtype=float32)\n"
     ]
    }
   ],
   "source": [
    "n = 10\n",
    "\n",
    "for i in range(n):\n",
    "    exec(\"x_{}=tf.constant(1.0)\".format(i+1), globals())\n",
    "    \n",
    "# for i in range(n):\n",
    "#     s = \"print('x_{0}: ', x_{0})\".format(i+1)\n",
    "#     exec(s, globals())\n",
    "    \n",
    "k = tf.constant([i+1.0 for i in range(n)])\n",
    "\n",
    "with tf.GradientTape(persistent=True) as g:\n",
    "    for i in range(n):\n",
    "        exec(\"g.watch(x_{})\".format(i+1), globals())\n",
    "    temp = tf.constant(0.)\n",
    "    for i in range(n):\n",
    "        # exec(\"temp = k[{}]*x_{}\".format(i, i+1), globals())\n",
    "        exec(\"temp += tf.multiply(k[{}], x_{})\".format(i, i+1), globals())\n",
    "        # print(temp)\n",
    "    z = tf.pow(temp,2)\n",
    "    # print(z)\n",
    "    \n",
    "\n",
    "for i in range(n):\n",
    "    s = \"dz_dx{0} = g.gradient(z, x_{0})\".format(i+1)\n",
    "    exec(s, globals())\n",
    "    s = \"print('dz_dx{0}: ', dz_dx{0})\".format(i+1)\n",
    "    exec(s, globals())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18b004c1",
   "metadata": {},
   "source": [
    "第二种方法：不批量命名变量，考虑为一个向量"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "fbfcc201",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tf.Tensor([ 110.  220.  330.  440.  550.  660.  770.  880.  990. 1100.], shape=(10,), dtype=float32)\n"
     ]
    }
   ],
   "source": [
    "n = 10  # 0, 1, 10,100,1000,10000---\n",
    "x = tf.ones((n))  # 注意参数是元组 (n)\n",
    "k = tf.constant([i+1.0 for i in range(n)])  # float\n",
    "with tf.GradientTape(persistent=True) as g:\n",
    "    g.watch(x)\n",
    "    z = tf.pow(tf.reduce_sum(tf.multiply(k, x)),2) # reduce_sum\n",
    "    \n",
    "dz_dx = g.gradient(z, x)\n",
    "print(dz_dx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "eedb642f",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(), dtype=float32, numpy=3025.0002>"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "z"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45a0c4d0",
   "metadata": {},
   "source": [
    "tf.reduce_sum 计算张量维度上元素的总和。如果忘记了，就错了。注意到这里z并不同，tensorflow会自动广播而不报错，需要小心谨慎。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "id": "13ea83e9",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tf.Tensor(\n",
      "[  2.        8.       18.       32.       50.       72.       98.00001\n",
      " 128.      162.      200.     ], shape=(10,), dtype=float32)\n"
     ]
    }
   ],
   "source": [
    "n = 10  # 0, 1, 10,100,1000,10000---\n",
    "x = tf.ones((n))  # 注意参数是元组 (n)\n",
    "k = tf.constant([i+1.0 for i in range(n)])  # float\n",
    "with tf.GradientTape(persistent=True) as g:\n",
    "    g.watch(x)\n",
    "    z = tf.pow(tf.multiply(k, x),2)\n",
    "    \n",
    "dz_dx = g.gradient(z, x)\n",
    "print(dz_dx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "id": "4a37bdd2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(10,), dtype=float32, numpy=\n",
       "array([ 1.      ,  4.      ,  9.      , 16.      , 24.999998, 36.      ,\n",
       "       48.999996, 64.      , 81.      , 99.99999 ], dtype=float32)>"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "id": "0d01161c",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tf.Tensor(\n",
      "[6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6.\n",
      " 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6.\n",
      " 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6.\n",
      " 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6. 6.\n",
      " 6. 6. 6. 6.], shape=(100,), dtype=float32) tf.Tensor([1200.], shape=(1,), dtype=float32)\n"
     ]
    }
   ],
   "source": [
    "n = 100\n",
    "x = tf.ones((n+1))  # 最后一维作为时间\n",
    "# x = tf.ones((n+n))\n",
    "(x,y) = x[:n],x[n:]\n",
    "\n",
    "with tf.GradientTape(persistent=True) as g:\n",
    "    g.watch(x)\n",
    "    g.watch(y)\n",
    "    z = (x+2*y)*(x+2*y)\n",
    "    \n",
    "dz_dx = g.gradient(z, x)\n",
    "dz_dy = g.gradient(z, y)\n",
    "print(dz_dx, dz_dy)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a4da349",
   "metadata": {},
   "source": [
    "矩阵、张量\n",
    "\n",
    "$\\large \\mathbf{z} = (\\mathbf{x}+2\\mathbf{y})^2 = \\mathbf{x}^2+4\\mathbf{y}^2+4\\mathbf{x}\\mathbf{y}$\n",
    "\n",
    "$\\large \\frac{\\partial \\mathbf{z}}{\\partial \\mathbf{x}}\\mid_{(\\mathbf{x},\\mathbf{y})=(\\mathbf{1},\\mathbf{0})}, \\frac{\\partial \\mathbf{z}}{\\partial \\mathbf{y}}\\mid_{(\\mathbf{x},\\mathbf{y})=(\\mathbf{1},\\mathbf{0})}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "id": "2589f3dd",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tf.Tensor(\n",
      "[[4. 4. 4. 4.]\n",
      " [4. 4. 4. 4.]\n",
      " [4. 4. 4. 4.]\n",
      " [4. 4. 4. 4.]], shape=(4, 4), dtype=float32) tf.Tensor(\n",
      "[[8. 8. 8. 8.]\n",
      " [8. 8. 8. 8.]\n",
      " [8. 8. 8. 8.]\n",
      " [8. 8. 8. 8.]], shape=(4, 4), dtype=float32)\n"
     ]
    }
   ],
   "source": [
    "x = tf.zeros((4,4))  # ;y = tf.ones((4,4))\n",
    "# x = tf.ones((4,4,4)) #;y = tf.ones((4,4,4))\n",
    "# y = tf.constant(0.0)\n",
    "y = tf.ones((4,4))\n",
    "\n",
    "with tf.GradientTape(persistent=True) as g:\n",
    "    g.watch(x)\n",
    "    g.watch(y)\n",
    "    z = (x+2*y)*(x+2*y)\n",
    "    \n",
    "dz_dx = g.gradient(z, x)\n",
    "dz_dy = g.gradient(z, y)\n",
    "print(dz_dx, dz_dy)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3575057",
   "metadata": {},
   "source": [
    "### n元函数的二阶偏导------黑塞矩阵"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a5f62de",
   "metadata": {},
   "source": [
    "**Hessian Matrix 黑塞矩阵**\n",
    "\n",
    "设 $Ω$ 是 $R^n$ 中的一个集合，$x=(x_1,x_2,\\dots,x_n)$ 表示 $Ω$ 上的点。假设 $u=u(x):Ω \\rightarrow R$ 是一个函数。对固定正整数 $k$，我们用符号 $D^k u$ 表示 $u$ 的所有 $k$ 阶偏导数 $\\large \\frac {\\partial^k u} {\\partial x_{i_1}\\partial x_{i_2} \\cdots \\partial x_{i_k}}$，其中 $(i_1,i_2,\\cdots,i_k)$ 是集合 $\\{1,2,\\cdots,n\\}$ 中 $k$ 个元素的任意排列。$D^k u$ 可以被看成是 $n^k$ 维欧氏空间 $R^{n^k}$ 上的向量。\n",
    "\n",
    "当 $k=1$ 时，我们称 $n$ 维向量 $\\large Du = (\\frac {\\partial u} {\\partial x_1}, \\frac {\\partial u} {\\partial x_2}, \\cdots, \\frac {\\partial u} {\\partial x_n})$ 为 $u$ 的**梯度**；\n",
    "\n",
    "当 $k=2$ 时，我们称 $n \\times n$ 矩阵\n",
    "$$\n",
    "\\Large\n",
    "D^2u = \\left[ \\begin{matrix}\n",
    "\\frac {\\partial^2 u} {\\partial x^2_1} &  \\frac {\\partial^2 u} {\\partial x_1 \\partial x_2} &  \\cdots &  \\frac {\\partial^2 u} {\\partial x_1 \\partial x_n} \\\\\n",
    "\\frac {\\partial^2 u} {\\partial x_2 \\partial x_1} &  \\frac {\\partial^2 u} {\\partial x^2_2} &  \\cdots &  \\frac {\\partial^2 u} {\\partial x_2 \\partial x_n} \\\\\n",
    "\\vdots & \\vdots & \\ddots & \\vdots                 \\\\\n",
    "\\frac {\\partial^2 u} {\\partial x_n \\partial x_1} &  \\frac {\\partial^2 u} {\\partial x_n \\partial x_2} &  \\cdots &  \\frac {\\partial^2 u} {\\partial x^2_n} \\\\\n",
    "\\end{matrix} \\right]\n",
    "$$\n",
    "为 $u$ 的**Hessian 矩阵**；"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0cb89dda",
   "metadata": {},
   "source": [
    "二元函数的二阶偏导\n",
    "\n",
    "$$\n",
    "\\large \\frac{\\partial^2}{\\partial x^2}(z = (x+2y)^2)\\mid_{(x,y)=(1,0)}, \\frac{\\partial^2}{\\partial y^2}(z = (x+2y)^2)\\mid_{(x,y)=(1,0)}, \\\\ \n",
    "\\large \\frac{\\partial^2}{\\partial x \\partial y}(z = (x+2y)^2)\\mid_{(x,y)=(1,0)}, \\frac{\\partial^2}{\\partial y \\partial x }(z = (x+2y)^2)\\mid_{(x,y)=(1,0)}.\\\\\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c492fea9",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tf.Tensor(2.0, shape=(), dtype=float32) tf.Tensor(4.0, shape=(), dtype=float32)\n",
      "tf.Tensor(2.0, shape=(), dtype=float32) tf.Tensor(4.0, shape=(), dtype=float32)\n",
      "tf.Tensor(4.0, shape=(), dtype=float32) tf.Tensor(8.0, shape=(), dtype=float32)\n"
     ]
    }
   ],
   "source": [
    "x = tf.constant(1.0)\n",
    "y = tf.constant(0.0)\n",
    "\n",
    "with tf.GradientTape(persistent=True) as g:\n",
    "    g.watch(x)\n",
    "    g.watch(y)\n",
    "    with tf.GradientTape(persistent=True) as gg:\n",
    "        gg.watch(x)\n",
    "        gg.watch(y)\n",
    "        z = (x+2*y)*(x+2*y)\n",
    "    dz_dx = gg.gradient(z, x)\n",
    "    dz_dy = gg.gradient(z, y)\n",
    "    \n",
    "d2z_dx2 = g.gradient(dz_dx, x)\n",
    "d2z_dxy = g.gradient(dz_dx, y)\n",
    "\n",
    "d2z_dy2 = g.gradient(dz_dy, y)\n",
    "d2z_dyx = g.gradient(dz_dy, x)\n",
    "\n",
    "print(dz_dx, dz_dy) # 一阶\n",
    "\n",
    "print(d2z_dx2, d2z_dxy) # 二阶\n",
    "print(d2z_dyx, d2z_dy2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "340d4c7a",
   "metadata": {},
   "source": [
    "三元函数的二阶偏导"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c6811249",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tf.Tensor(12.0, shape=(), dtype=float32) tf.Tensor(24.0, shape=(), dtype=float32) tf.Tensor(36.0, shape=(), dtype=float32)\n",
      "tf.Tensor(2.0, shape=(), dtype=float32) tf.Tensor(4.0, shape=(), dtype=float32) tf.Tensor(6.0, shape=(), dtype=float32)\n",
      "tf.Tensor(4.0, shape=(), dtype=float32) tf.Tensor(8.0, shape=(), dtype=float32) tf.Tensor(12.0, shape=(), dtype=float32)\n",
      "tf.Tensor(6.0, shape=(), dtype=float32) tf.Tensor(12.0, shape=(), dtype=float32) tf.Tensor(18.0, shape=(), dtype=float32)\n"
     ]
    }
   ],
   "source": [
    "x = tf.constant(1.0)\n",
    "y = tf.constant(1.0)\n",
    "m = tf.constant(1.0)\n",
    "with tf.GradientTape(persistent=True) as g:\n",
    "    g.watch(x)\n",
    "    g.watch(y)\n",
    "    g.watch(m)\n",
    "    with tf.GradientTape(persistent=True) as gg:\n",
    "        gg.watch(x)\n",
    "        gg.watch(y)\n",
    "        gg.watch(m)\n",
    "        z = (x+2*y+3*m)*(x+2*y+3*m)\n",
    "    dz_dx = gg.gradient(z, x)\n",
    "    dz_dy = gg.gradient(z, y)\n",
    "    dz_dm = gg.gradient(z, m)\n",
    "    \n",
    "d2z_dx2 = g.gradient(dz_dx, x)\n",
    "d2z_dxy = g.gradient(dz_dx, y)\n",
    "d2z_dxm = g.gradient(dz_dx, m)\n",
    "\n",
    "d2z_dy2 = g.gradient(dz_dy, y)\n",
    "d2z_dyx = g.gradient(dz_dy, x)\n",
    "d2z_dym = g.gradient(dz_dy, m)\n",
    "\n",
    "d2z_dm2 = g.gradient(dz_dm, m)\n",
    "d2z_dmx = g.gradient(dz_dm, x)\n",
    "d2z_dmy = g.gradient(dz_dm, y)\n",
    "\n",
    "\n",
    "print(dz_dx, dz_dy, dz_dm)\n",
    "print(d2z_dx2, d2z_dxy, d2z_dxm)\n",
    "print(d2z_dyx, d2z_dy2, d2z_dym)\n",
    "print(d2z_dmx, d2z_dmy, d2z_dm2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "982982d7",
   "metadata": {},
   "source": [
    "n元函数的二阶偏导"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "16171ecf",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(<tf.Tensor: shape=(5,), dtype=float32, numpy=array([ 30.,  60.,  90., 120., 150.], dtype=float32)>,\n",
       " <tf.Tensor: shape=(5, 5), dtype=float32, numpy=\n",
       " array([[ 2.,  4.,  6.,  8., 10.],\n",
       "        [ 4.,  8., 12., 16., 20.],\n",
       "        [ 6., 12., 18., 24., 30.],\n",
       "        [ 8., 16., 24., 32., 40.],\n",
       "        [10., 20., 30., 40., 50.]], dtype=float32)>)"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = 5\n",
    "\n",
    "for i in range(n):\n",
    "    exec(\"x_{}=tf.constant(1.0)\".format(i+1), globals())\n",
    "# for i in range(n):\n",
    "#     s = \"print('x_{0}: ', x_{0})\".format(i+1)\n",
    "#     exec(s, globals())\n",
    "\n",
    "k = tf.constant([i+1.0 for i in range(n)])\n",
    "\n",
    "with tf.GradientTape(persistent=True) as g:\n",
    "    for i in range(n):\n",
    "        exec(\"g.watch(x_{})\".format(i+1), globals())\n",
    "        \n",
    "    with tf.GradientTape(persistent=True) as gg:\n",
    "        for i in range(n):\n",
    "            exec(\"gg.watch(x_{})\".format(i+1), globals())\n",
    "            \n",
    "        temp = tf.constant(0.)\n",
    "        for i in range(n):\n",
    "            exec(\"temp += tf.multiply(k[{}], x_{})\".format(i, i+1), globals())\n",
    "            # print(temp)\n",
    "        z = tf.pow(temp,2)\n",
    "        # print(z)\n",
    "        \n",
    "    gradient_vector = np.zeros((n))\n",
    "    for i in range(n):\n",
    "        s = \"dz_dx{0} = gg.gradient(z, x_{0})\".format(i+1)\n",
    "        exec(s, globals())\n",
    "        s = \"gradient_vector[{0}] = dz_dx{1}\".format(i,i+1)\n",
    "        exec(s, globals())\n",
    "        # s = \"print('dz_dx{0}: ', dz_dx{0})\".format(i+1)\n",
    "        # exec(s, globals())\n",
    "    gradient_vector = tf.constant(gradient_vector, dtype=tf.float32)\n",
    "        \n",
    "hessian_matrix = np.zeros((n,n))\n",
    "for i in range(n):\n",
    "    for j in range(n):\n",
    "        s = \"d2z_dx{0}x{1} = g.gradient(dz_dx{0}, x_{1})\".format(i+1,j+1)\n",
    "        exec(s, globals())\n",
    "        s = \"hessian_matrix[{0}][{1}] = d2z_dx{2}x{3}\".format(i,j,i+1,j+1)\n",
    "        exec(s, globals())\n",
    "        \n",
    "        # s = \"print('d2z_dx{0}x{1}: ', d2z_dx{0}x{1}.numpy(),end='\\t')\".format(i+1,j+1)\n",
    "        # exec(s, globals())\n",
    "    # print()\n",
    "    \n",
    "hessian_matrix = tf.constant(hessian_matrix, dtype=tf.float32)\n",
    "\n",
    "gradient_vector, hessian_matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c92c4d0d",
   "metadata": {},
   "source": [
    "d2z_dx1xn这些变量不定义会返回 nan，但是定义了又没有什么用吧？\n",
    "\n",
    "```python\n",
    "n = 5\n",
    "\n",
    "for i in range(n):\n",
    "    exec(\"x_{}=tf.constant(1.0)\".format(i+1), globals())\n",
    "\n",
    "k = tf.constant([i+1.0 for i in range(n)])\n",
    "\n",
    "with tf.GradientTape(persistent=True) as g:\n",
    "    for i in range(n):\n",
    "        exec(\"g.watch(x_{})\".format(i+1), globals())\n",
    "        \n",
    "    with tf.GradientTape(persistent=True) as gg:\n",
    "        for i in range(n):\n",
    "            exec(\"gg.watch(x_{})\".format(i+1), globals())\n",
    "            \n",
    "        temp = tf.constant(0.)\n",
    "        for i in range(n):\n",
    "            exec(\"temp += tf.multiply(k[{}], x_{})\".format(i, i+1), globals())\n",
    "        z = tf.pow(temp,2)\n",
    "        \n",
    "    gradient_vector = np.zeros((n))\n",
    "    for i in range(n):\n",
    "        s = \"gradient_vector[{0}] = gg.gradient(z, x_{1})\".format(i,i+1)\n",
    "        exec(s, globals())\n",
    "    gradient_vector = tf.constant(gradient_vector, dtype=tf.float32)\n",
    "        \n",
    "hessian_matrix = np.zeros((n,n))\n",
    "for i in range(n):\n",
    "    for j in range(n):\n",
    "        s = \"hessian_matrix[{0}][{1}] = g.gradient(gradient_vector[{2}], x_{3})\".format(i,j,i,j+1)\n",
    "        exec(s, globals())\n",
    "hessian_matrix = tf.constant(hessian_matrix, dtype=tf.float32)\n",
    "\n",
    "gradient_vector, hessian_matrix\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6ae3f16",
   "metadata": {},
   "source": [
    "**GradientTape.jacobian**\n",
    "\n",
    "tf 得到标量函数的二次导数,但这种模式并不能通用于生成黑塞矩阵, 因为 `GradientTape.gradient` 只计算标量的梯度。要构造黑塞矩阵，使用 `GradientTape.jacobian` 方法\n",
    "\n",
    "对于未连接的梯度，得到 0 而不是 None 会比较方便。您可以使用 unconnected_gradients 参数来决定具有未连接的梯度时返回的内容\n",
    "\n",
    "tape.gradient(z, x, unconnected_gradients=tf.UnconnectedGradients.ZERO)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "id": "9d7e3dbe",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(<tf.Tensor: shape=(5,), dtype=float32, numpy=array([ 30.,  60.,  90., 120., 150.], dtype=float32)>,\n",
       " <tf.Tensor: shape=(5, 5), dtype=float32, numpy=\n",
       " array([[ 2.,  4.,  6.,  8., 10.],\n",
       "        [ 4.,  8., 12., 16., 20.],\n",
       "        [ 6., 12., 18., 24., 30.],\n",
       "        [ 8., 16., 24., 32., 40.],\n",
       "        [10., 20., 30., 40., 50.]], dtype=float32)>)"
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = 5\n",
    "x = tf.Variable(tf.ones((n)))  # 注意参数是元组 (n)\n",
    "k = tf.constant([i+1.0 for i in range(n)])  # float\n",
    "\n",
    "with tf.GradientTape(persistent=True) as g:\n",
    "    with tf.GradientTape(persistent=True) as gg:\n",
    "        z = tf.pow(tf.reduce_sum(tf.multiply(k, x)),2)\n",
    "    dz_dx = gg.jacobian(z, x)\n",
    "\n",
    "# d2z_dx2 = g.gradient(dz_dx, x)  # 只能得到标量函数的二次导数\n",
    "d2z_dx2 = g.jacobian(dz_dx, x)    # 生成黑塞矩阵\n",
    "\n",
    "dz_dx, d2z_dx2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "id": "1a1258e0",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(<tf.Tensor: shape=(5,), dtype=float32, numpy=array([ 675., 1350., 2025., 2700., 3375.], dtype=float32)>,\n",
       " <tf.Tensor: shape=(5, 5), dtype=float32, numpy=\n",
       " array([[  90.,  180.,  270.,  360.,  450.],\n",
       "        [ 180.,  360.,  540.,  720.,  900.],\n",
       "        [ 270.,  540.,  810., 1080., 1350.],\n",
       "        [ 360.,  720., 1080., 1440., 1800.],\n",
       "        [ 450.,  900., 1350., 1800., 2250.]], dtype=float32)>,\n",
       " <tf.Tensor: shape=(5, 5, 5), dtype=float32, numpy=\n",
       " array([[[  6.,  12.,  18.,  24.,  30.],\n",
       "         [ 12.,  24.,  36.,  48.,  60.],\n",
       "         [ 18.,  36.,  54.,  72.,  90.],\n",
       "         [ 24.,  48.,  72.,  96., 120.],\n",
       "         [ 30.,  60.,  90., 120., 150.]],\n",
       " \n",
       "        [[ 12.,  24.,  36.,  48.,  60.],\n",
       "         [ 24.,  48.,  72.,  96., 120.],\n",
       "         [ 36.,  72., 108., 144., 180.],\n",
       "         [ 48.,  96., 144., 192., 240.],\n",
       "         [ 60., 120., 180., 240., 300.]],\n",
       " \n",
       "        [[ 18.,  36.,  54.,  72.,  90.],\n",
       "         [ 36.,  72., 108., 144., 180.],\n",
       "         [ 54., 108., 162., 216., 270.],\n",
       "         [ 72., 144., 216., 288., 360.],\n",
       "         [ 90., 180., 270., 360., 450.]],\n",
       " \n",
       "        [[ 24.,  48.,  72.,  96., 120.],\n",
       "         [ 48.,  96., 144., 192., 240.],\n",
       "         [ 72., 144., 216., 288., 360.],\n",
       "         [ 96., 192., 288., 384., 480.],\n",
       "         [120., 240., 360., 480., 600.]],\n",
       " \n",
       "        [[ 30.,  60.,  90., 120., 150.],\n",
       "         [ 60., 120., 180., 240., 300.],\n",
       "         [ 90., 180., 270., 360., 450.],\n",
       "         [120., 240., 360., 480., 600.],\n",
       "         [150., 300., 450., 600., 750.]]], dtype=float32)>)"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = 5\n",
    "x = tf.Variable(tf.ones((n)))  # 注意参数是元组 (n)\n",
    "k = tf.constant([i+1.0 for i in range(n)])  # float\n",
    "with tf.GradientTape(persistent=True) as g:\n",
    "    with tf.GradientTape(persistent=True) as gg:\n",
    "        with tf.GradientTape(persistent=True) as ggg:\n",
    "            z = tf.pow(tf.reduce_sum(tf.multiply(k, x)),3)\n",
    "        dz_dx = ggg.jacobian(z, x)\n",
    "    \n",
    "    # d2z_dx2 = g.gradient(dz_dx, x)  # 只能得到标量函数的二次导数\n",
    "    d2z_dx2 = gg.jacobian(dz_dx, x)    # 生成黑塞矩阵\n",
    "\n",
    "d3z_dx3 = g.jacobian(d2z_dx2, x)    # 生成黑塞矩阵\n",
    "\n",
    "dz_dx, d2z_dx2, d3z_dx3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "id": "f194031e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "id": "4815139d",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(<tf.Tensor: shape=(2, 2), dtype=float32, numpy=\n",
       " array([[12., 24.],\n",
       "        [12., 24.]], dtype=float32)>,\n",
       " <tf.Tensor: shape=(2, 2, 2, 2), dtype=float32, numpy=\n",
       " array([[[[2., 4.],\n",
       "          [2., 4.]],\n",
       " \n",
       "         [[4., 8.],\n",
       "          [4., 8.]]],\n",
       " \n",
       " \n",
       "        [[[2., 4.],\n",
       "          [2., 4.]],\n",
       " \n",
       "         [[4., 8.],\n",
       "          [4., 8.]]]], dtype=float32)>)"
      ]
     },
     "execution_count": 77,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = 2  # 0, 1, 10,100,1000,10000---\n",
    "m = 2\n",
    "x = tf.ones((m, n))  # 注意参数是元组 (n)\n",
    "k = tf.constant([i+1.0 for i in range(n)])  # float\n",
    "\n",
    "with tf.GradientTape(persistent=True) as g:\n",
    "    g.watch(x)\n",
    "    with tf.GradientTape(persistent=True) as gg:\n",
    "        gg.watch(x)\n",
    "        z = tf.pow(tf.reduce_sum(tf.multiply(k, x)),2)\n",
    "    dz_dx = gg.gradient(z, x)\n",
    "\n",
    "d2z_dx2 = g.jacobian(dz_dx, x)\n",
    "\n",
    "dz_dx, d2z_dx2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bfb21a84",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "f289ffb4",
   "metadata": {},
   "source": [
    "### m元函数的n阶偏导、n阶全微分\n",
    "\n",
    "*假设函数在定义域上无穷次可微*。\n",
    "\n",
    "对一个m元函数$\\large y = f(x_1,x_2,\\dots,x_m)$求偏导得到一阶偏导数$df$，函数在定义域内关于$x_i$的偏导数记作$\\large \\frac{\\partial f}{\\partial x_i} 或 f_{x_i}$, 函数在点$x_0$关于$x_i$的偏导数$\\large \\frac{\\partial f}{\\partial x_i} (x_0) = \\frac{\\partial f}{\\partial x_i}\\mid_{x_0}$，函数在点$x_0$的方向导数(关于$x_1,x_2,\\dots,x_m$的偏导数)是**梯度**，记作$$\\large grad f(x_0) = \\nabla f(x_0) = (\\frac{\\partial f}{\\partial x_1}(x_0), \\frac{\\partial f}{\\partial x_2}(x_0),\\dots, \\frac{\\partial f}{\\partial x_m}(x_0)).$$\n",
    "\n",
    "函数$f$在点$x_0$的**全微分**记作$df(x_0)$,即\n",
    "$$\\large df(x_0) = \\sum^m_{i=1} \\frac{\\partial f}{\\partial x_i}(x_0)dx_i = \\sum^m_{i=1} f_{x_i}(x_0)dx_i .$$\n",
    "\n",
    "对一个m元函数$\\large y = f(x_1,x_2,\\dots,x_m)$关于$x_i$求偏导得到一阶偏导数$\\large \\frac{\\partial f}{\\partial x_i}$，对一阶偏导数再求偏导数得到二阶偏导数$\\large \\frac{\\partial^2 f}{\\partial x_j \\partial x_i}$，对二阶偏导数再求偏导数得到三阶偏导数$\\large \\frac{\\partial^3 f}{\\partial x_k \\partial x_j \\partial x_i }$，可得n阶偏导数为$\\large \\frac{\\partial^n f}{\\partial x_{i_n} \\dots \\partial x_{i_2} \\partial x_{i_1} } = f_{x_{i_n}\\dots x_{i_2}x_{i_1}} .$\n",
    "\n",
    "m元函数$f$在点$x_0$的一阶全微分记作$df(x_0)$,即$\\large df(x_0) = \\sum^m_{i=1} \\frac{\\partial f}{\\partial x_i}(x_0)dx_i = \\sum^m_{i=1} f_{x_i}(x_0)dx_i .$\n",
    "\n",
    "m元函数$f$在点$x_0$的二阶全微分记作$d^2 f(x_0)$,是一个以对称矩阵$\\large (f_{x_i x_j}(x_0))_{m\\times m}$为系数矩阵的二次型，叫黑塞矩阵，即\n",
    "$$\n",
    "\\large \\boldsymbol{H} f(x_0) = (f_{x_i x_j}(x_0))_{m\\times m} = \n",
    "\\left[ \\begin{matrix}\n",
    "\\frac {\\partial^2 f} {\\partial x^2_1} (x_0) &  \\frac {\\partial^2 f} {\\partial x_1 \\partial x_2} (x_0) &  \\cdots &  \\frac {\\partial^2 f} {\\partial x_1 \\partial x_m} (x_0) \\\\\n",
    "\\frac {\\partial^2 f} {\\partial x_2 \\partial x_1} (x_0) &  \\frac {\\partial^2 f} {\\partial x^2_2} (x_0)&  \\cdots &  \\frac {\\partial^2 f} {\\partial x_2 \\partial x_m} (x_0)\\\\\n",
    "\\vdots & \\vdots & \\ddots & \\vdots                 \\\\\n",
    "\\frac {\\partial^2 f} {\\partial x_m \\partial x_1} (x_0)&  \\frac {\\partial^2 f} {\\partial x_m \\partial x_2} (x_0) &  \\cdots &  \\frac {\\partial^2 f} {\\partial x^2_m} (x_0) \\\\\n",
    "\\end{matrix} \\right].\n",
    "$$\n",
    "\n",
    "m元函数$f$在点$x_0$的n阶全微分记作$d^n f(x_0)$，它的系数构成一个n阶张量$T_nf(x_0)$，它的形状是$(m,m,\\dots,m)$，一共有n个m(在同构的意义下，第零阶张量为标量，第一阶张量为向量， 第二阶张量为矩阵)。经典的方法把张量视为多维数组，张量的元素是数组中的值。\n",
    "$$\\large \\boldsymbol{T_n} f(x_0) = (f_{x_{i_n}\\dots x_{i_2}x_{i_1}} (x))_{(m \\times m \\times \\dots \\times m)_n} .$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5037ecf4",
   "metadata": {},
   "source": [
    "### m元复合函数的偏导数\n",
    "\n",
    "一个n维m元向量函数\n",
    "$$\\large \\boldsymbol{F}(\\boldsymbol{x}) = (F_1(x_1,x_2,\\dots,x_m), F_2(x_1,x_2,\\dots,x_m), \\dots, F_n(x_1,x_2,\\dots,x_m))$$\n",
    "一个可微函数$g$, 考虑m元复合函数记为 $\\large g \\circ \\boldsymbol{F} = g(\\boldsymbol{F}(x))$。\n",
    "\n",
    "\n",
    "向量函数$\\boldsymbol{F}(x)$的**雅可比矩阵**，记作\n",
    "$$\n",
    "\\large\n",
    "\\frac{D \\boldsymbol{F}(x)}{D x} = \\left[ \\begin{matrix}\n",
    "\\frac {\\partial F_1(x)} {\\partial x_1} &  \\frac {\\partial F_1(x)} {\\partial x_2} &  \\cdots &  \\frac {\\partial F_1(x)} {\\partial x_m} \\\\\n",
    "\\frac {\\partial F_2(x)} {\\partial x_1} &  \\frac {\\partial F_2(x)} {\\partial x_2} &  \\cdots &  \\frac {\\partial F_2(x)} {\\partial x_m} \\\\\n",
    "\\vdots & \\vdots & \\ddots & \\vdots                 \\\\\n",
    "\\frac {\\partial F_n(x)} {\\partial x_1} &  \\frac {\\partial F_n(x)} {\\partial x_2} &  \\cdots &  \\frac {\\partial F_n(x)} {\\partial x_m} \\\\\n",
    "\\end{matrix} \\right].\n",
    "$$\n",
    "\n",
    "则复合函数 $g \\circ \\boldsymbol{F}$ 的偏导数为：$$\\nabla (g(\\boldsymbol{F}(x))) = \\nabla g(y) \\mit_{y=\\boldsymbol{F}(x)} \\frac{D \\boldsymbol{F}(x)}{D x}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64838b96",
   "metadata": {},
   "source": [
    "## 张量求导\n",
    "\n",
    "\n",
    "#### 对标量求导\n",
    "\n",
    "1. 标量对标量求导，结果是个标量\n",
    "事实上就是一元函数求一阶导$\\frac{dy}{dx}$，通过一元函数的n阶导求n+1阶导。\n",
    "2. 向量对标量求导，结果是个向量\n",
    "事实上就是向量的每一个元素对标量求导，n个一元函数(线性无关的n元函数)求一阶导$\\frac{d y_i}{d x_i}$。\n",
    "3. 矩阵对标量求导，结果是个矩阵\n",
    "事实上也就是矩阵的每一个元素对标量求导，$m \\times n$个一元函数求一阶导$\\frac{d y_{ij}}{d x_{ij}}$。\n",
    "4. n阶张量对标量求导，结果是个矩阵\n",
    "事实上也就是n阶张量的每一个元素对标量求导，$m_1 \\times m_2 \\times \\dots \\times m_n$个一元函数求一阶导$\\frac{d y_{i_{1}...i_{n}}}{d x_{i_{1}...i_{n}}}$。\n",
    "\n",
    "####  对向量求导\n",
    "\n",
    "1. 标量对向量求导，结果是向量\n",
    "事实上这就是所谓的Gradient，记为n元(标量)函数的梯度 $\\nabla f$ ；也可以看成多元函数求一阶偏导数。\n",
    "2. 向量对向量求导，结果是矩阵\n",
    "这个当然也是gradient，当然这准确的说应该叫matrix gradient. 可以看成对多元函数一阶偏导数求二阶偏导数，又如向量函数$\\boldsymbol{F}(x)$的雅可比矩阵。\n",
    "3. 矩阵对向量求导，结果是个三阶张量\n",
    "对多元函数二阶偏导数求三阶偏导数。\n",
    "4. n阶张量对向量求导，结果是n+1阶张量\n",
    "对多元函数n阶偏导数求n+1阶偏导数。\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "####  对矩阵求导\n",
    "\n",
    "1. 标量对矩阵求导，结果是矩阵\n",
    "事实上这一类，主要是考虑一类标量函数对矩阵的导数，一般是det,trace,log(det)等等\n",
    "2. 向量对矩阵求导，结果是3阶张量\n",
    "对多元函数n阶偏导数求n+1阶偏导数。\n",
    "3. 矩阵对矩阵求导，结果是4阶张量\n",
    "对多元函数n阶偏导数求n+1阶偏导数。\n",
    "4. n阶张量对矩阵求导，结果是n+2阶张量\n",
    "对多元函数n阶偏导数求n+1阶偏导数。\n",
    "\n",
    "#### 对张量求导\n",
    "\n",
    "**结论：n阶张量对m阶张量求导，结果是n+m阶张量。**\n",
    "\n",
    "标量是0阶张量，向量是1阶张量，矩阵是2阶张量。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40ad1e70",
   "metadata": {},
   "source": [
    "### tf张量求导\n",
    "\n",
    "GradientTape.jacobian, GradientTape.batch_jacobian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "id": "f043c585",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(2, 2, 2, 2), dtype=float32, numpy=\n",
       "array([[[[2., 0.],\n",
       "         [0., 0.]],\n",
       "\n",
       "        [[0., 4.],\n",
       "         [0., 0.]]],\n",
       "\n",
       "\n",
       "       [[[0., 0.],\n",
       "         [6., 0.]],\n",
       "\n",
       "        [[0., 0.],\n",
       "         [0., 8.]]]], dtype=float32)>"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "with tf.GradientTape() as g:\n",
    "    x = tf.constant([[1., 2.], [3., 4.]], dtype=tf.float32)\n",
    "    g.watch(x)\n",
    "    y = x * x\n",
    "jacobian = g.jacobian(y, x)\n",
    "jacobian # [[[2,  0], [0,  4]], [[6,  0], [0,  8]]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "id": "7c599ddd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(2, 2, 2), dtype=float32, numpy=\n",
       "array([[[2., 0.],\n",
       "        [0., 4.]],\n",
       "\n",
       "       [[6., 0.],\n",
       "        [0., 8.]]], dtype=float32)>"
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "with tf.GradientTape() as g:\n",
    "    x = tf.constant([[1., 2.], [3., 4.]], dtype=tf.float32)\n",
    "    g.watch(x)\n",
    "    y = x * x\n",
    "batch_jacobian = g.batch_jacobian(y, x)\n",
    "batch_jacobian # [[[2,  0], [0,  4]], [[6,  0], [0,  8]]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "id": "5801c46e",
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(<tf.Tensor: shape=(5,), dtype=float32, numpy=array([ 675., 1350., 2025., 2700., 3375.], dtype=float32)>,\n",
       " <tf.Tensor: shape=(5, 5), dtype=float32, numpy=\n",
       " array([[  90.,  180.,  270.,  360.,  450.],\n",
       "        [ 180.,  360.,  540.,  720.,  900.],\n",
       "        [ 270.,  540.,  810., 1080., 1350.],\n",
       "        [ 360.,  720., 1080., 1440., 1800.],\n",
       "        [ 450.,  900., 1350., 1800., 2250.]], dtype=float32)>,\n",
       " <tf.Tensor: shape=(5, 5, 5), dtype=float32, numpy=\n",
       " array([[[  6.,  12.,  18.,  24.,  30.],\n",
       "         [ 12.,  24.,  36.,  48.,  60.],\n",
       "         [ 18.,  36.,  54.,  72.,  90.],\n",
       "         [ 24.,  48.,  72.,  96., 120.],\n",
       "         [ 30.,  60.,  90., 120., 150.]],\n",
       " \n",
       "        [[ 12.,  24.,  36.,  48.,  60.],\n",
       "         [ 24.,  48.,  72.,  96., 120.],\n",
       "         [ 36.,  72., 108., 144., 180.],\n",
       "         [ 48.,  96., 144., 192., 240.],\n",
       "         [ 60., 120., 180., 240., 300.]],\n",
       " \n",
       "        [[ 18.,  36.,  54.,  72.,  90.],\n",
       "         [ 36.,  72., 108., 144., 180.],\n",
       "         [ 54., 108., 162., 216., 270.],\n",
       "         [ 72., 144., 216., 288., 360.],\n",
       "         [ 90., 180., 270., 360., 450.]],\n",
       " \n",
       "        [[ 24.,  48.,  72.,  96., 120.],\n",
       "         [ 48.,  96., 144., 192., 240.],\n",
       "         [ 72., 144., 216., 288., 360.],\n",
       "         [ 96., 192., 288., 384., 480.],\n",
       "         [120., 240., 360., 480., 600.]],\n",
       " \n",
       "        [[ 30.,  60.,  90., 120., 150.],\n",
       "         [ 60., 120., 180., 240., 300.],\n",
       "         [ 90., 180., 270., 360., 450.],\n",
       "         [120., 240., 360., 480., 600.],\n",
       "         [150., 300., 450., 600., 750.]]], dtype=float32)>)"
      ]
     },
     "execution_count": 79,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = 5\n",
    "x = tf.Variable(tf.ones((n)))  # 注意参数是元组 (n)\n",
    "k = tf.constant([i+1.0 for i in range(n)])  # float\n",
    "with tf.GradientTape(persistent=True) as g:\n",
    "    with tf.GradientTape(persistent=True) as gg:\n",
    "        with tf.GradientTape(persistent=True) as ggg:\n",
    "            z = tf.pow(tf.reduce_sum(tf.multiply(k, x)),3)\n",
    "        dz_dx = ggg.jacobian(z, x)\n",
    "    \n",
    "    # d2z_dx2 = g.gradient(dz_dx, x)  # 只能得到标量函数的二次导数\n",
    "    d2z_dx2 = gg.jacobian(dz_dx, x)    # 生成黑塞矩阵\n",
    "\n",
    "d3z_dx3 = g.jacobian(d2z_dx2, x)    # 生成三阶张量\n",
    "\n",
    "dz_dx, d2z_dx2, d3z_dx3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8f302d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def diff(func:str, x, n:int, ifexec=False):\n",
    "    if n<1:\n",
    "        global y\n",
    "        exec(func, globals())\n",
    "    else:\n",
    "        with tf.GradientTape(persistent=True) as g:\n",
    "            g.watch(x)\n",
    "            diff(func, x, n-1,ifexec)\n",
    "        try:\n",
    "            y = g.jacobian(y, x)\n",
    "        except TypeError:\n",
    "            pass\n",
    "    # print(n, y)\n",
    "    if y== None: y = tf.constant(0.)\n",
    "    if ifexec:\n",
    "        # exec executes a Python statement stored in a string or file\n",
    "        exec(\"d{0}y_dx{0} = y\".format(n), globals())\n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "155eb48c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'x_1': 1, 'x_2': 2, 'x_3': 3}"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = [1,2,3]\n",
    "dic = {\"x_1\": n[0],\"x_2\": n[1],\"x_3\": n[2],}\n",
    "\n",
    "dic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "id": "5201370c",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "temp = tf.constant(0.)\n",
      "for i in range(n):\n",
      "    exec(\"temp += tf.multiply(k[{}], x_{})\".format(i, i+1), globals())\n",
      "    # print(temp)\n",
      "z = tf.pow(temp,2)\n",
      " False\n",
      "x_1 1\n",
      "x_2 0\n",
      "x_3 1\n",
      "3 2\n"
     ]
    }
   ],
   "source": [
    "def diff2(func:str, kwargs:dict, ifexec=False):\n",
    "    print(func, ifexec)\n",
    "    order_number = 0\n",
    "    n_dim =  0\n",
    "    for key, value in kwargs.items():\n",
    "        print(key, value)\n",
    "        order_number += value\n",
    "        n_dim += 1\n",
    "    print(n_dim, order_number)\n",
    "    \n",
    "    \n",
    "\n",
    "n = [1,0,1]\n",
    "\n",
    "for i in range(len(n)):\n",
    "    exec(\"x_{}=tf.constant(1.0)\".format(i+1), globals())\n",
    "\n",
    "func = \"\"\"\n",
    "temp = tf.constant(0.)\n",
    "for i in range(n):\n",
    "    exec(\"temp += tf.multiply(k[{}], x_{})\".format(i, i+1), globals())\n",
    "    # print(temp)\n",
    "z = tf.pow(temp,2)\n",
    "\"\"\"\n",
    "dic = {\"x_1\": n[0],\"x_2\": n[1],\"x_3\": n[2],}\n",
    "diff2(func=func, kwargs=dic)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "ca6e9e20",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = [1,2,3]\n",
    "\n",
    "len(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9208ff6b",
   "metadata": {},
   "outputs": [],
   "source": [
    "diff(func,{x_1,n[0]}, {x_2,n[0]},{x_3,n[0]},)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac529d20",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 5\n",
    "\n",
    "for i in range(n):\n",
    "    exec(\"x_{}=tf.constant(1.0)\".format(i+1), globals())\n",
    "# for i in range(n):\n",
    "#     s = \"print('x_{0}: ', x_{0})\".format(i+1)\n",
    "#     exec(s, globals())\n",
    "\n",
    "k = tf.constant([i+1.0 for i in range(n)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7768767",
   "metadata": {},
   "outputs": [],
   "source": [
    "diff(func, x, nlist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "id": "e68e9263",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6\n",
      "6\n"
     ]
    }
   ],
   "source": [
    "func = \"\"\"\n",
    "a = 1\n",
    "b = 2\n",
    "exec(\"c = 3\")\n",
    "exec(\"print(a+b+c)\")\n",
    "print(a+b+c)\n",
    "\"\"\"\n",
    "\n",
    "exec(func)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "id": "fec43743",
   "metadata": {},
   "outputs": [],
   "source": [
    "def diff(func:str, x, n:int, ifexec=False):\n",
    "    \"\"\"\n",
    "    func: a string, such as func=\"y = x * x\";\n",
    "    x: tf.FloatTensor, such as x=tf.constant(1.0);\n",
    "    n: a int number, Order of derivative, such as n=1;\n",
    "    ifexec: \n",
    "    return dny_dxn: the value of order of derivative.\n",
    "    example: \n",
    "        x = tf.constant(1.0)\n",
    "        n = 3\n",
    "        func = \"y = x * x * x\"\n",
    "        diff(func, x, n)\n",
    "    \"\"\"\n",
    "    if n<1:\n",
    "        global y\n",
    "        exec(func, globals())\n",
    "    else:\n",
    "        with tf.GradientTape(persistent=True) as g:\n",
    "            g.watch(x)\n",
    "            diff(func, x, n-1,ifexec)\n",
    "        try:\n",
    "            y = g.gradient(y, x)\n",
    "        except TypeError:\n",
    "            pass\n",
    "    # print(n, y)\n",
    "    if y== None: y = tf.constant(0.)\n",
    "    if ifexec:\n",
    "        # exec executes a Python statement stored in a string or file\n",
    "        exec(\"d{0}y_dx{0} = y\".format(n), globals())\n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "id": "d202c3cb",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(4,), dtype=float32, numpy=array([0.31705797, 0.18140507, 1.71776   , 3.513605  ], dtype=float32)>"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = tf.constant([1.,2.,3.,4.])\n",
    "n = 2\n",
    "func = \"y = 2.*tf.sin(x)+tf.pow(x,2)\"\n",
    "diff(func=func, ifexec=False, x=x,n=n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "4e91bcc2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "y = 2*tf.sin(x) False\n",
      "x 5\n"
     ]
    }
   ],
   "source": [
    "def diff2(func:str, kwargs, ifexec=False, ):\n",
    "    print(func, ifexec)\n",
    "    for key, value in kwargs.items():\n",
    "        print(key, value)\n",
    "        \n",
    "x = tf.constant(0.0)\n",
    "n = 5\n",
    "func = \"y = 2*tf.sin(x)\"\n",
    "diff2(func=func, kwargs={\"x\": n})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63859f46",
   "metadata": {},
   "source": [
    "https://reference.wolfram.com/language/ref/D.html.zh\n",
    "\n",
    "D[f, {x, n}, {y, n}]\n",
    "\n",
    "D[Sin[x]^10, {x, 4}]\n",
    "\n",
    "D[Sin[x y]/(x^2 + y^2), x, y]\n",
    "\n",
    "D[x^2 E^(5 y), {x, 2}, {y, 3}]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "id": "4b9f61c8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: shape=(), dtype=float32, numpy=2.0>"
      ]
     },
     "execution_count": 97,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = tf.constant(0.0)\n",
    "n = 5\n",
    "func = \"y = 2*tf.sin(x)\"\n",
    "diff(func, x, n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b82450d",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = tf.constant(0.0)\n",
    "n = 5\n",
    "func = \"y = 2*tf.sin(x)\"\n",
    "diff(func, x, n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65bd93c2",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "9f2bee74",
   "metadata": {},
   "source": [
    "## SymPy\n",
    "\n",
    "SymPy 是一个符号计算的Python库.\n",
    "\n",
    "\n",
    "```\n",
    "conda install sympy\n",
    "pip install sympy\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "id": "47d71aa2",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import symbols, diff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "id": "aa8be285",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x**2 + x*y + y**2 + 2\n",
      "9\n",
      "2*x + y\n",
      "4\n",
      "x + 2*y\n",
      "5\n"
     ]
    }
   ],
   "source": [
    "x, y = symbols('x, y')\n",
    "\n",
    "z = x**2+y**2+x*y+2\n",
    "print(z)\n",
    "result = z.subs({x: 1, y: 2})   # 用数值分别对x、y进行替换\n",
    "print(result)\n",
    "\n",
    "dx = diff(z, x)   # 对x求偏导\n",
    "print(dx)\n",
    "result = dx.subs({x: 1, y: 2})\n",
    "print(result)\n",
    "\n",
    "dy = diff(z, y)   # 对y求偏导\n",
    "print(dy)\n",
    "result = dy.subs({x: 1, y: 2})\n",
    "print(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39dabad8",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c819af07",
   "metadata": {},
   "outputs": [],
   "source": [
    "geom = dde.geometry.Interval(-1, 1)\n",
    "timedomain = dde.geometry.TimeDomain(0, 1)\n",
    "geomtime = dde.geometry.GeometryXTime(geom, timedomain)\n",
    "\n",
    "lambda_1 = 1\n",
    "lambda_2 = 0.0025\n",
    "d = 100\n",
    "\n",
    "def pde(x, y):\n",
    "    dy_x = tf.gradients(y, x)[0]\n",
    "    dy_x, dy_t = dy_x[:, 0:d], dy_x[:, d:]\n",
    "    dy_xx = tf.gradients(dy_x, x)[0][:, 0:d]\n",
    "    dy_xxx = tf.gradients(dy_xx, x)[0][:, 0:d]\n",
    "    return dy_t + lambda_1*y*dy_x +  lambda_2*dy_xxx\n",
    "\n",
    "bc = dde.icbc.DirichletBC(geomtime, lambda x: 0, lambda _, on_boundary: on_boundary)\n",
    "ic = dde.icbc.IC(geomtime, lambda x: np.cos(np.pi * x[:, 0:d]), lambda _, on_initial: on_initial)\n",
    "\n",
    "data = dde.data.TimePDE(geomtime, pde, [bc, ic], num_domain=8000, num_boundary=400, num_initial=800)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7e35cfc",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"Backend supported: tensorflow.compat.v1, tensorflow, pytorch\"\"\"\n",
    "import deepxde as dde\n",
    "import numpy as np\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "id": "dda43d23",
   "metadata": {},
   "outputs": [],
   "source": [
    "geom = dde.geometry.Interval(-1, 1)\n",
    "timedomain = dde.geometry.TimeDomain(0, 1)\n",
    "geomtime = dde.geometry.GeometryXTime(geom, timedomain)\n",
    "\n",
    "lambda_1 = 1\n",
    "lambda_2 = 0.0025\n",
    "d = 100\n",
    "\n",
    "def pde(x, y):\n",
    "    dy_x = tf.gradients(y, x)[0]\n",
    "    dy_x, dy_t = dy_x[:, 0:d], dy_x[:, d:]\n",
    "    dy_xx = tf.gradients(dy_x, x)[0][:, 0:d]\n",
    "    dy_xxx = tf.gradients(dy_xx, x)[0][:, 0:d]\n",
    "    return dy_t + lambda_1*y*dy_x +  lambda_2*dy_xxx\n",
    "\n",
    "bc = dde.icbc.DirichletBC(geomtime, lambda x: 0, lambda _, on_boundary: on_boundary)\n",
    "ic = dde.icbc.IC(geomtime, lambda x: np.cos(np.pi * x[:, 0:d]), lambda _, on_initial: on_initial)\n",
    "\n",
    "data = dde.data.TimePDE(geomtime, pde, [bc, ic], num_domain=8000, num_boundary=400, num_initial=800)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 155,
   "id": "9b406fae",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10400, 2)"
      ]
     },
     "execution_count": 155,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.train_next_batch(64)[0].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "id": "c215833d",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Compiling model...\n",
      "Building feed-forward neural network...\n",
      "'build' took 0.053891 s\n",
      "\n"
     ]
    },
    {
     "ename": "RuntimeError",
     "evalue": "IC func should return an array of shape N by 1 for a single component.Use argument 'component' for different components.",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mRuntimeError\u001b[0m                              Traceback (most recent call last)",
      "Input \u001b[1;32mIn [142]\u001b[0m, in \u001b[0;36m<cell line: 4>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      1\u001b[0m net \u001b[38;5;241m=\u001b[39m dde\u001b[38;5;241m.\u001b[39mnn\u001b[38;5;241m.\u001b[39mFNN([d\u001b[38;5;241m+\u001b[39m\u001b[38;5;241m1\u001b[39m] \u001b[38;5;241m+\u001b[39m [\u001b[38;5;241m20\u001b[39m] \u001b[38;5;241m*\u001b[39m \u001b[38;5;241m3\u001b[39m \u001b[38;5;241m+\u001b[39m [\u001b[38;5;241m1\u001b[39m], \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mtanh\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mGlorot normal\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[0;32m      2\u001b[0m model \u001b[38;5;241m=\u001b[39m dde\u001b[38;5;241m.\u001b[39mModel(data, net)\n\u001b[1;32m----> 4\u001b[0m \u001b[43mmodel\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mcompile\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43madam\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mlr\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m1e-3\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[0;32m      5\u001b[0m model\u001b[38;5;241m.\u001b[39mtrain(epochs\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m15000\u001b[39m)\n\u001b[0;32m      6\u001b[0m model\u001b[38;5;241m.\u001b[39mcompile(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mL-BFGS\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n",
      "File \u001b[1;32mG:\\Anaconda3\\envs\\py3.8\\lib\\site-packages\\deepxde\\utils\\internal.py:22\u001b[0m, in \u001b[0;36mtiming.<locals>.wrapper\u001b[1;34m(*args, **kwargs)\u001b[0m\n\u001b[0;32m     19\u001b[0m \u001b[38;5;129m@wraps\u001b[39m(f)\n\u001b[0;32m     20\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mwrapper\u001b[39m(\u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs):\n\u001b[0;32m     21\u001b[0m     ts \u001b[38;5;241m=\u001b[39m timeit\u001b[38;5;241m.\u001b[39mdefault_timer()\n\u001b[1;32m---> 22\u001b[0m     result \u001b[38;5;241m=\u001b[39m \u001b[43mf\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m     23\u001b[0m     te \u001b[38;5;241m=\u001b[39m timeit\u001b[38;5;241m.\u001b[39mdefault_timer()\n\u001b[0;32m     24\u001b[0m     \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;132;01m%r\u001b[39;00m\u001b[38;5;124m took \u001b[39m\u001b[38;5;132;01m%f\u001b[39;00m\u001b[38;5;124m s\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;241m%\u001b[39m (f\u001b[38;5;241m.\u001b[39m\u001b[38;5;18m__name__\u001b[39m, te \u001b[38;5;241m-\u001b[39m ts))\n",
      "File \u001b[1;32mG:\\Anaconda3\\envs\\py3.8\\lib\\site-packages\\deepxde\\model.py:108\u001b[0m, in \u001b[0;36mModel.compile\u001b[1;34m(self, optimizer, lr, loss, metrics, decay, loss_weights, external_trainable_variables)\u001b[0m\n\u001b[0;32m    105\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mexternal_trainable_variables \u001b[38;5;241m=\u001b[39m external_trainable_variables\n\u001b[0;32m    107\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m backend_name \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mtensorflow.compat.v1\u001b[39m\u001b[38;5;124m\"\u001b[39m:\n\u001b[1;32m--> 108\u001b[0m     \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_compile_tensorflow_compat_v1\u001b[49m\u001b[43m(\u001b[49m\u001b[43mlr\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mloss_fn\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdecay\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mloss_weights\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m    109\u001b[0m \u001b[38;5;28;01melif\u001b[39;00m backend_name \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mtensorflow\u001b[39m\u001b[38;5;124m\"\u001b[39m:\n\u001b[0;32m    110\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_compile_tensorflow(lr, loss_fn, decay, loss_weights)\n",
      "File \u001b[1;32mG:\\Anaconda3\\envs\\py3.8\\lib\\site-packages\\deepxde\\model.py:130\u001b[0m, in \u001b[0;36mModel._compile_tensorflow_compat_v1\u001b[1;34m(self, lr, loss_fn, decay, loss_weights)\u001b[0m\n\u001b[0;32m    127\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39msaver \u001b[38;5;241m=\u001b[39m tf\u001b[38;5;241m.\u001b[39mtrain\u001b[38;5;241m.\u001b[39mSaver(max_to_keep\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m)\n\u001b[0;32m    129\u001b[0m \u001b[38;5;66;03m# Data losses\u001b[39;00m\n\u001b[1;32m--> 130\u001b[0m losses \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mdata\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mlosses\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mnet\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mtargets\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mnet\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43moutputs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mloss_fn\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[0;32m    131\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(losses, \u001b[38;5;28mlist\u001b[39m):\n\u001b[0;32m    132\u001b[0m     losses \u001b[38;5;241m=\u001b[39m [losses]\n",
      "File \u001b[1;32mG:\\Anaconda3\\envs\\py3.8\\lib\\site-packages\\deepxde\\data\\pde.py:154\u001b[0m, in \u001b[0;36mPDE.losses\u001b[1;34m(self, targets, outputs, loss, model)\u001b[0m\n\u001b[0;32m    152\u001b[0m     beg, end \u001b[38;5;241m=\u001b[39m bcs_start[i], bcs_start[i \u001b[38;5;241m+\u001b[39m \u001b[38;5;241m1\u001b[39m]\n\u001b[0;32m    153\u001b[0m     \u001b[38;5;66;03m# The same BC points are used for training and testing.\u001b[39;00m\n\u001b[1;32m--> 154\u001b[0m     error \u001b[38;5;241m=\u001b[39m \u001b[43mbc\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43merror\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mtrain_x\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmodel\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mnet\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43minputs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43moutputs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mbeg\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mend\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m    155\u001b[0m     losses\u001b[38;5;241m.\u001b[39mappend(loss[\u001b[38;5;28mlen\u001b[39m(error_f) \u001b[38;5;241m+\u001b[39m i](bkd\u001b[38;5;241m.\u001b[39mzeros_like(error), error))\n\u001b[0;32m    156\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m losses\n",
      "File \u001b[1;32mG:\\Anaconda3\\envs\\py3.8\\lib\\site-packages\\deepxde\\icbc\\initial_conditions.py:32\u001b[0m, in \u001b[0;36mIC.error\u001b[1;34m(self, X, inputs, outputs, beg, end)\u001b[0m\n\u001b[0;32m     30\u001b[0m values \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mfunc(X, beg, end)\n\u001b[0;32m     31\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m bkd\u001b[38;5;241m.\u001b[39mndim(values) \u001b[38;5;241m>\u001b[39m \u001b[38;5;241m0\u001b[39m \u001b[38;5;129;01mand\u001b[39;00m bkd\u001b[38;5;241m.\u001b[39mshape(values)[\u001b[38;5;241m1\u001b[39m] \u001b[38;5;241m!=\u001b[39m \u001b[38;5;241m1\u001b[39m:\n\u001b[1;32m---> 32\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mRuntimeError\u001b[39;00m(\n\u001b[0;32m     33\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mIC func should return an array of shape N by 1 for a single component.\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m     34\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mUse argument \u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mcomponent\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m for different components.\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m     35\u001b[0m     )\n\u001b[0;32m     36\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m outputs[beg:end, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mcomponent : \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mcomponent \u001b[38;5;241m+\u001b[39m \u001b[38;5;241m1\u001b[39m] \u001b[38;5;241m-\u001b[39m values\n",
      "\u001b[1;31mRuntimeError\u001b[0m: IC func should return an array of shape N by 1 for a single component.Use argument 'component' for different components."
     ]
    }
   ],
   "source": [
    "\n",
    "\n",
    "\n",
    "net = dde.nn.FNN([d+1] + [20] * 3 + [1], \"tanh\", \"Glorot normal\")\n",
    "model = dde.Model(data, net)\n",
    "\n",
    "model.compile(\"adam\", lr=1e-3)\n",
    "model.train(epochs=15000)\n",
    "model.compile(\"L-BFGS\")\n",
    "losshistory, train_state = model.train()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f90db538",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Using backend: tensorflow.compat.v1\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING:tensorflow:From G:\\Anaconda3\\envs\\py3.8\\lib\\site-packages\\tensorflow\\python\\compat\\v2_compat.py:101: disable_resource_variables (from tensorflow.python.ops.variable_scope) is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "non-resource variables are not supported in the long term\n",
      "WARNING:tensorflow:From G:\\Anaconda3\\envs\\py3.8\\lib\\site-packages\\deepxde\\nn\\initializers.py:118: The name tf.keras.initializers.he_normal is deprecated. Please use tf.compat.v1.keras.initializers.he_normal instead.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import deepxde as dde\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "956f555f",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "G:\\Anaconda3\\envs\\py3.8\\lib\\site-packages\\skopt\\sampler\\sobol.py:246: UserWarning: The balance properties of Sobol' points require n to be a power of 2. 0 points have been previously generated, then: n=0+2542=2542. \n",
      "  warnings.warn(\"The balance properties of Sobol' points require \"\n",
      "G:\\Anaconda3\\envs\\py3.8\\lib\\site-packages\\skopt\\sampler\\sobol.py:246: UserWarning: The balance properties of Sobol' points require n to be a power of 2. 0 points have been previously generated, then: n=0+82=82. \n",
      "  warnings.warn(\"The balance properties of Sobol' points require \"\n",
      "G:\\Anaconda3\\envs\\py3.8\\lib\\site-packages\\skopt\\sampler\\sobol.py:246: UserWarning: The balance properties of Sobol' points require n to be a power of 2. 0 points have been previously generated, then: n=0+162=162. \n",
      "  warnings.warn(\"The balance properties of Sobol' points require \"\n"
     ]
    }
   ],
   "source": [
    "import deepxde as dde\n",
    "import numpy as np\n",
    "\n",
    "def pde(x, y):\n",
    "    dy_x = dde.grad.jacobian(y, x, i=0, j=0)\n",
    "    dy_t = dde.grad.jacobian(y, x, i=0, j=1)\n",
    "    dy_xx = dde.grad.hessian(y, x, i=0, j=0)\n",
    "    return dy_t + y * dy_x - 0.01 / np.pi * dy_xx\n",
    "\n",
    "geom = dde.geometry.Interval(-1, 1)\n",
    "timedomain = dde.geometry.TimeDomain(0, 0.99)\n",
    "geomtime = dde.geometry.GeometryXTime(geom, timedomain)\n",
    "\n",
    "bc = dde.icbc.DirichletBC(geomtime, lambda x: 0, lambda _, on_boundary: on_boundary)\n",
    "ic = dde.icbc.IC(\n",
    "    geomtime, lambda x: -np.sin(np.pi * x[:, 0:1]), lambda _, on_initial: on_initial\n",
    ")\n",
    "\n",
    "data = dde.data.TimePDE(\n",
    "    geomtime, pde, [bc, ic], num_domain=2540, num_boundary=80, num_initial=160\n",
    ")\n",
    "\n",
    "for i in [1, 10, 128]:\n",
    "    print(data.train_next_batch(batch_size=i)[0].shape)\n",
    "    print(data.test(batch_size=i)[0].shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "491d2fac",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(9200, 2)\n",
      "(9200, 2)\n",
      "(9200, 2)\n",
      "(9200, 2)\n",
      "(9200, 2)\n",
      "(9200, 2)\n"
     ]
    }
   ],
   "source": [
    "for i in [1, 10, 128]:\n",
    "    print(data.train_next_batch(batch_size=i)[0].shape)\n",
    "    print(data.test(batch_size=i)[0].shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "fd96ac0a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(9200, 2)\n",
      "(9200, 2)\n"
     ]
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4ec549bd",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3020, 2)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "3600c365",
   "metadata": {},
   "outputs": [],
   "source": [
    "import deepxde as dde\n",
    "import numpy as np\n",
    "from scipy.io import loadmat\n",
    "# Import tf if using backend tensorflow.compat.v1 or tensorflow\n",
    "from deepxde.backend import tf\n",
    "# Import torch if using backend pytorch\n",
    "# import torch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "29b618e4",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "G:\\Anaconda3\\envs\\py3.8\\lib\\site-packages\\skopt\\sampler\\sobol.py:246: UserWarning: The balance properties of Sobol' points require n to be a power of 2. 0 points have been previously generated, then: n=0+8002=8002. \n",
      "  warnings.warn(\"The balance properties of Sobol' points require \"\n",
      "G:\\Anaconda3\\envs\\py3.8\\lib\\site-packages\\skopt\\sampler\\sobol.py:246: UserWarning: The balance properties of Sobol' points require n to be a power of 2. 0 points have been previously generated, then: n=0+402=402. \n",
      "  warnings.warn(\"The balance properties of Sobol' points require \"\n",
      "G:\\Anaconda3\\envs\\py3.8\\lib\\site-packages\\skopt\\sampler\\sobol.py:246: UserWarning: The balance properties of Sobol' points require n to be a power of 2. 0 points have been previously generated, then: n=0+802=802. \n",
      "  warnings.warn(\"The balance properties of Sobol' points require \"\n"
     ]
    }
   ],
   "source": [
    "\n",
    "geom = dde.geometry.Interval(-1, 1)\n",
    "timedomain = dde.geometry.TimeDomain(0, 1)\n",
    "geomtime = dde.geometry.GeometryXTime(geom, timedomain)\n",
    "d = 0.001\n",
    "\n",
    "def pde(x, y):\n",
    "    dy_t = dde.grad.jacobian(y, x, i=0, j=1)\n",
    "    dy_xx = dde.grad.hessian(y, x, i=0, j=0)\n",
    "    return dy_t - d * dy_xx - 5 * (y - y**3)\n",
    "\n",
    "# Hard restraints on initial + boundary conditions\n",
    "# Backend tensorflow.compat.v1 or tensorflow\n",
    "def output_transform(x, y):\n",
    "    return x[:, 0:1]**2 * tf.cos(np.pi * x[:, 0:1]) + x[:, 1:2] * (1 - x[:, 0:1]**2) * y\n",
    "\n",
    "# Backend pytorch\n",
    "# def output_transform(x, y):\n",
    "#     return x[:, 0:1]**2 * torch.cos(np.pi * x[:, 0:1]) + x[:, 1:2] * (1 - x[:, 0:1]**2) * y\n",
    "\n",
    "data = dde.data.TimePDE(geomtime, pde, [], num_domain=8000, num_boundary=400, num_initial=800)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "4e2d8a2d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(9200, 2)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.train_next_batch(batch_size=7)[0].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "2b90b94f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(9200, 2)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.test()[0].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "83b39c48",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Using backend: tensorflow.compat.v1\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING:tensorflow:From G:\\Anaconda3\\envs\\py3.8\\lib\\site-packages\\tensorflow\\python\\compat\\v2_compat.py:101: disable_resource_variables (from tensorflow.python.ops.variable_scope) is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "non-resource variables are not supported in the long term\n",
      "WARNING:tensorflow:From G:\\Anaconda3\\envs\\py3.8\\lib\\site-packages\\deepxde\\nn\\initializers.py:118: The name tf.keras.initializers.he_normal is deprecated. Please use tf.compat.v1.keras.initializers.he_normal instead.\n",
      "\n",
      "(3020, 2)\n",
      "(3020, 2)\n",
      "(3020, 2)\n",
      "(3020, 2)\n",
      "(3020, 2)\n",
      "(3020, 2)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "G:\\Anaconda3\\envs\\py3.8\\lib\\site-packages\\skopt\\sampler\\sobol.py:246: UserWarning: The balance properties of Sobol' points require n to be a power of 2. 0 points have been previously generated, then: n=0+2542=2542. \n",
      "  warnings.warn(\"The balance properties of Sobol' points require \"\n",
      "G:\\Anaconda3\\envs\\py3.8\\lib\\site-packages\\skopt\\sampler\\sobol.py:246: UserWarning: The balance properties of Sobol' points require n to be a power of 2. 0 points have been previously generated, then: n=0+82=82. \n",
      "  warnings.warn(\"The balance properties of Sobol' points require \"\n",
      "G:\\Anaconda3\\envs\\py3.8\\lib\\site-packages\\skopt\\sampler\\sobol.py:246: UserWarning: The balance properties of Sobol' points require n to be a power of 2. 0 points have been previously generated, then: n=0+162=162. \n",
      "  warnings.warn(\"The balance properties of Sobol' points require \"\n"
     ]
    }
   ],
   "source": [
    "import deepxde as dde\n",
    "import numpy as np\n",
    "\n",
    "def pde(x, y):\n",
    "    dy_x = dde.grad.jacobian(y, x, i=0, j=0)\n",
    "    dy_t = dde.grad.jacobian(y, x, i=0, j=1)\n",
    "    dy_xx = dde.grad.hessian(y, x, i=0, j=0)\n",
    "    return dy_t + y * dy_x - 0.01 / np.pi * dy_xx\n",
    "\n",
    "geom = dde.geometry.Interval(-1, 1)\n",
    "timedomain = dde.geometry.TimeDomain(0, 0.99)\n",
    "geomtime = dde.geometry.GeometryXTime(geom, timedomain)\n",
    "\n",
    "bc = dde.icbc.DirichletBC(geomtime, lambda x: 0, lambda _, on_boundary: on_boundary)\n",
    "ic = dde.icbc.IC(\n",
    "    geomtime, lambda x: -np.sin(np.pi * x[:, 0:1]), lambda _, on_initial: on_initial\n",
    ")\n",
    "\n",
    "data = dde.data.TimePDE(\n",
    "    geomtime, pde, [bc, ic], num_domain=2540, num_boundary=80, num_initial=160\n",
    ")\n",
    "\n",
    "for i in [1, 10, 128]:\n",
    "    print(data.train_next_batch(batch_size=i)[0].shape)\n",
    "    print(data.test(batch_size=i)[0].shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f1118ee",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a9799b11",
   "metadata": {},
   "outputs": [],
   "source": [
    "# https://github.com/lululxvi/deepxde/blob/master/examples/pinn_forward/fractional_Poisson_3d.py\n",
    "\"\"\"Backend supported: tensorflow.compat.v1\"\"\"\n",
    "import deepxde as dde\n",
    "import numpy as np\n",
    "from deepxde.backend import tf\n",
    "from scipy.special import gamma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7f38ebc6",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "G:\\Anaconda3\\envs\\py3.8\\lib\\site-packages\\skopt\\sampler\\sobol.py:246: UserWarning: The balance properties of Sobol' points require n to be a power of 2. 0 points have been previously generated, then: n=0+258=258. \n",
      "  warnings.warn(\"The balance properties of Sobol' points require \"\n",
      "G:\\Anaconda3\\envs\\py3.8\\lib\\site-packages\\skopt\\sampler\\sobol.py:246: UserWarning: The balance properties of Sobol' points require n to be a power of 2. 0 points have been previously generated, then: n=0+3=3. \n",
      "  warnings.warn(\"The balance properties of Sobol' points require \"\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Warning: mesh step size 0.010000 is larger than the boundary distance 0.001304.\n"
     ]
    }
   ],
   "source": [
    "alpha = 1.8\n",
    "\n",
    "\n",
    "def fpde(x, y, int_mat):\n",
    "    \"\"\"\\int_theta D_theta^alpha u(x)\"\"\"\n",
    "    if isinstance(int_mat, (list, tuple)) and len(int_mat) == 3:\n",
    "        int_mat = tf.SparseTensor(*int_mat)\n",
    "        lhs = tf.sparse_tensor_dense_matmul(int_mat, y)\n",
    "    else:\n",
    "        lhs = tf.matmul(int_mat, y)\n",
    "    lhs = lhs[:, 0]\n",
    "    lhs *= gamma((1 - alpha) / 2) * gamma((3 + alpha) / 2) / (2 * np.pi ** 2)\n",
    "    x = x[: tf.size(lhs)]\n",
    "    rhs = (\n",
    "        2 ** alpha\n",
    "        * gamma(2 + alpha / 2)\n",
    "        * gamma((3 + alpha) / 2)\n",
    "        / gamma(3 / 2)\n",
    "        * (1 - (1 + alpha / 3) * tf.reduce_sum(x ** 2, axis=1))\n",
    "    )\n",
    "    return lhs - rhs\n",
    "\n",
    "\n",
    "def func(x):\n",
    "    return (np.abs(1 - np.linalg.norm(x, axis=1, keepdims=True) ** 2)) ** (\n",
    "        1 + alpha / 2\n",
    "    )\n",
    "\n",
    "\n",
    "geom = dde.geometry.Sphere([0, 0, 0], 1)\n",
    "bc = dde.icbc.DirichletBC(geom, func, lambda _, on_boundary: on_boundary)\n",
    "\n",
    "data = dde.data.FPDE(\n",
    "    geom,\n",
    "    fpde,\n",
    "    alpha,\n",
    "    bc,\n",
    "    [8, 8, 100],\n",
    "    num_domain=256,\n",
    "    num_boundary=1,\n",
    "    solution=func,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "73437563",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1251989, 3)\n",
      "(1251988, 3)\n",
      "(1251989, 3)\n",
      "(1251988, 3)\n",
      "(1251989, 3)\n",
      "(1251988, 3)\n"
     ]
    }
   ],
   "source": [
    "for i in [1, 10, 128]:\n",
    "    print(data.train_next_batch(batch_size=i)[0].shape)\n",
    "    print(data.test(batch_size=i)[0].shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9fb6be9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "net = dde.nn.FNN([3] + [20] * 4 + [1], \"tanh\", \"Glorot normal\")\n",
    "net.apply_output_transform(\n",
    "    lambda x, y: (1 - tf.reduce_sum(x ** 2, axis=1, keepdims=True)) * y\n",
    ")\n",
    "\n",
    "model = dde.Model(data, net)\n",
    "model.compile(\"adam\", lr=1e-3)\n",
    "losshistory, train_state = model.train(epochs=10000)\n",
    "dde.saveplot(losshistory, train_state, issave=False, isplot=True)\n",
    "\n",
    "X = geom.random_points(10000)\n",
    "y_true = func(X)\n",
    "y_pred = model.predict(X)\n",
    "print(\"L2 relative error:\", dde.metrics.l2_relative_error(y_true, y_pred))\n",
    "np.savetxt(\"test.dat\", np.hstack((X, y_true, y_pred)))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py3.8",
   "language": "python",
   "name": "py3.8"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
